Object-Based Greenhouse Horticultural Crop Identification from Multi-Temporal Satellite Imagery : A Case Study in Almeria , Spain

Greenhouse detection and mapping via remote sensing is a complex task, which has already been addressed in numerous studies. In this research, the innovative goal relies on the identification of greenhouse horticultural crops that were growing under plastic coverings on 30 September 2013. To this end, object-based image analysis (OBIA) and a decision tree classifier (DT) were applied to a set consisting of eight Landsat 8 OLI images collected from May to November 2013. Moreover, a single WorldView-2 satellite image acquired on 30 September 2013, was also used as a data source. In this approach, basic spectral information, textural features and several vegetation indices (VIs) derived from Landsat 8 and WorldView-2 multi-temporal satellite data were computed on previously segmented image objects in order to identify four of the most popular autumn crops cultivated under greenhouse in Almería, Spain (i.e., tomato, pepper, cucumber and aubergine). The best classification accuracy (81.3% overall accuracy) was achieved by using the full set of Landsat 8 time series. These results were considered good in the case of tomato and pepper crops, being significantly worse for cucumber and aubergine. These OPEN ACCESS Remote Sens. 2015, 7 7379 results were hardly improved by adding the information of the WorldView-2 image. The most important information for correct classification of different crops under greenhouses was related to the greenhouse management practices and not the spectral properties of the crops themselves.


Introduction
The practice of using plastic materials to advance the first harvest and increase the yield of horticultural crops has been steadily increasing during the past 60 years in the entire world [1].Perhaps nowhere this process is more visible than around Almería, southern Spain, where agriculture under plastic has become a key economic driver.However, this productive sector is being affected by an unstable and changing geopolitical situation of the markets.In this regard, facts, such as the entry into force of the new agricultural agreement between the European Union (EU) and The Kingdom of Morocco adopted in 2012, which liberalizes EU-Morocco trade in agriculture (Council Decision of 8 March 2012, 2012/497/EU), or, more recently, the Russian veto of horticultural products from the EU, are causing important changes in the market prices of products grown under plastic-covered greenhouses, which significantly affect farmers and agribusiness.In order to alleviate these changes in market prices, each agricultural cooperative in Almería is carrying out a crop acreage plan to be applied by their farmers.However, a more globalized planning of horticultural production would be desirable.Certainly, the possibility of knowing the crops that are being grown under greenhouses along an agricultural campaign, both focused on our local productive sector and our direct competitors, would help decision making in order to keep prices up, thus maximizing returns to farmers.In that sense, useful information for decision makers would include cropland extent, crop types, crop rotation dynamics, and so forth [2].This sort of information over space and time can be repeatedly, consistently and accurately acquired by methods relying on remote sensing without exception [3].In fact, remote sensing can significantly contribute to providing a timely and accurate picture of the agricultural sector, as it is very suitable for gathering information over large areas with high revisit frequency [4].
According to several authors [5][6][7][8][9][10], greenhouse mapping via remote sensing is a challenging task due to the fact that the spectral signature of the plastic-covered greenhouse can change drastically.In this way, different plastic materials with varying thickness, transparency, ultraviolet and infrared reflection and transmission properties, additives, age and colors are used in greenhouse coverings.Moreover, plastic sheets may be painted white (whitewashed) during summer to protect plants from excessive radiation and to reduce the heat inside the greenhouse.Finally, as plastic sheets are semi-transparent, the changing reflectance of the crops underneath them affects the greenhouse spectral signal reaching the sensor.
The identification of outdoor crops by using Landsat time series has been already addressed by many authors [22,23].Recently, an innovative methodology, which combines multi-temporal moderate resolution remote sensing data, the OBIA approach and the decision tree (DT) classifier algorithm, has been proposed [24,25].It should be noted that not only spatial and spectral resolutions of remotely-sensed data used in this type of investigation are very important, but also the temporal resolution of data is crucial.In this sense, Hao et al. [26] tested the potential of time series merged from multi-source satellite imagery for outdoor crop classification.
This paper presents the first attempt to identify greenhouse horticultural crops growing under plastic coverings by applying OBIA and decision tree techniques to multi-temporal remote sensing data.So far, and to the best knowledge of the authors, this is the first research work dealing with the identification of protected cropping.To this end, dense Landsat 8 OLI time series and a single WorldView-2 satellite image have been used as data sources.More specific subjects regarding the greenhouse crop classification method addressed in this work are: (1) to study the influence on the classification results of both the temporal and geometric resolution of satellite imagery; (2) to evaluate the most important and useful features for both crop and greenhouse covering (basic spectral information, textural features and several VIs); and (3) to detect the optimal temporal windows for the satellite imagery time series, as well as crucial time periods.

Study Area
This investigation was conducted in Almería, southern Spain, which has become the site of the greatest concentration of greenhouses in the world, known as the "Sea of Plastic" or the "Poniente" region.The study area comprises a rectangle area of about 8000 ha centered on the WGS84 geographic coordinates of 36.7824°N and 2.6867°W (Figure 1).
In the 2013-2014 agricultural year (from July 2013 to June 2014), the intensive agriculture in Almería was mainly dedicated to tomato (26% of total area), pepper (22%), zucchini (16%), cucumber (11%), aubergine (4.5%) and green bean (3%).The rest of the area was cultivated with other crops, such as melon, watermelon or Chinese cabbage.It is important to note that two crop cycles are usually grown in a greenhouse during an agricultural season in Almería with the greenhouses' empty time being little.The "first crop" in Almería is often planted between June and September and the second one between December and February.Therefore, the different combinations and possibilities of the cropping calendar in a greenhouse are extremely high.It is worth noting that this work will focus on the identification of crops grown in a greenhouse during the first crop or summer planting corresponding to the 2013-2014 agricultural season.
Moreover, the types of greenhouse plastic covers in the study area are very variable.The most common materials are polyethylene films (e.g., low density, long life, thermic, with or without additives) and ethylene-vinyl acetate copolymer, also known as EVA.Furthermore, both materials can present different thickness (180 μm or 200 μm) and colors (white, yellow or green).
A single WV2 image taken on 30 September 2013, over the study area was acquired.It was collected in Ortho Ready Standard Level-2A (ORS2A) format, containing both PAN and MS imagery.This satellite image had a mean off-nadir view angle of 11.8°, a mean collection azimuth of 38.2° and 0% of cloud cover.The final product GSD was of 0.4 m and 1.6 m for PAN and MS images, respectively.WV2 ORS2A format presents both radiometric and geometric corrections.These types of images are georeferenced to a cartographic projection using a surface with constant height and include the corresponding RPC sensor camera model and metadata file.The delivered products were ordered with a dynamic range of 11-bit and without the application of the dynamic range adjustment preprocessing (i.e., maintaining absolute radiometric accuracy and full dynamic range for scientific applications).
From this WV2 ORS2A bundle image, a pan-sharpened image with 0.4-m GSD was attained by means of the PANSHARP module included in Geomatica v. 2014 (PCI Geomatics, Richmond Hill, ON, Canada).The coordinates of seven ground control points (GCPs) and 32 independent check points (ICPs) were obtained by the differential global positioning system (DGPS) using a total GPS Topcon HiPer PRO station working in real-time kinematic mode (RTK).The ground points (both GCPs and ICPs) were measured with reference to the European Terrestrial Reference System 1989 (ETRS89) and UTM Zone 30 projection.A pan-sharpened RGB orthoimage with 0.4-m GSD was generated by using the seven GCPs to compute the sensor model based on rational functions refined by a zero order transformation in the image space (RPC0).A medium resolution 10-m grid spacing DEM with a vertical accuracy of 1.34 m (root mean square error; RMSE), provided by the Andalusia Government, was used to carry out the orthorectification process.The planimetric accuracy (RMSE2D) attained on the orthorectified pan-sharpened image and measured on 32 ICPs was 0.59 m.
Furthermore, an MS orthoimage with 1.6-m GSD and containing the full 8-band spectral information was also produced.The same seven GCPs, RPC0 model and DEM were used to attain the MS orthoimage.It is worth noting that atmospheric correction is especially important in those cases where multi-temporal or multi-sensor images are analyzed [24,25,[27][28][29][30].Thus, in this case and as a prior step to the orthorectification process, the original WV2 MS image was atmospherically corrected by using the ATCOR (atmospheric correction) module included in Geomatica v. 2014.This absolute atmospheric correction algorithm involves the conversion of the original raw digital numbers to ground reflectance values, and it is based on the MODTRAN (MODerate resolution atmospheric TRANsmission) radiative transfer code [31].Finally, an atmospherically-corrected WV2 MS orthoimage with a RMSE2D of 2.20 m was generated.
A dense time series of the Landsat 8's OLI sensor was downloaded, totally free and via the Internet, from the U.S. Landsat archive [32].It was composed of eight free of cloud multi-temporal Landsat 8 scenes type Level 1 Terrain Corrected (L1T) in the footprint of Path 200 and Row 34, each one comprising the whole working area.The eight Landsat 8 scenes were taken with a dynamic range of 12 bits from May to November in 2013, and they were dated as follow: 23 May, 8 June, 24 June, 10 July, 26 July, 11 August, 12 September and 15 November.
In order to use the aforementioned time series of Landsat 8 for greenhouse horticultural crop identification, it was necessary to carry out some previous processing steps.The first one was to clip the original images, both PAN and MS, according to the study area.It is noteworthy that spatial resolutions equal to or better than 8-16 m are recommended by Levin et al. [6] to work on greenhouse detection via remote sensing.Thereby, and after clipping, eight pan-sharpened images were attained by fusing the information of PAN and MS Landsat 8 bands through applying the PANSHARP module included in Geomatica v. 2014.The pan-sharpened images were generated with a GSD of 15 m, 12-bit depth and seven bands (C, B, G, R, NIR, SWIR1 and SWIR2).
The next processing step consisted of performing atmospheric correction by applying the ATCOR module (Geomatica v. 2014) to the eight Landsat 8 pan-sharpened images.The goal was to attain ground reflectance values.At this point, it was checked that ground reflectance values attained over greenhouses on the corrected pan-sharpened images were not significantly different from those achieved on the atmospherically corrected Landsat MS images.The results showed that the coefficient of determination (R 2 ) between ground reflectance values of pan-sharpened and MS images was always higher than 0.99 for each of the seven studied bands.
The last process involved the rectification and co-registration of the eight images.Although all acquired Landsat 8 L1T images have been previously orthorectified and corrected for terrain relief [33], a subsequent co-registration is needed, since an extremely accurate spatial matching among multi-temporal images is crucial in this type of work [34,35] for both pixel and object-based image analysis.The pan-sharpened RGB WV2 orthoimage with a GSD of 0.4 m was used as the reference image.Each pan-sharpened Landsat 8 image was rectified with this reference image using 50 ground control points evenly distributed over the working area and a first order polynomial transformation.The resulting RMSE2D for the geometrically-corrected images ranged from 4.42 m to 6.95 m.

Reference Greenhouses
The reference objects or ground truth for this work consisted of 694 individual greenhouses located in the study area and managed by Vicasol, an agricultural cooperative of Almería.Vicasol provided valuable information related to those crops growing under plastic coverings corresponding to these reference greenhouses in autumn 2013, such as type of crop, date of planting and date of removal (Table 1).We can see that these data fit well with the typical transplant pattern in the study area (Table 2).This information consisted of an Excel file, where the cadastral reference of each greenhouse was provided.It is important to note that the Vicasol data only vaguely reflect the distribution of the main crops for the whole of Almería.This fact should be taken into account in further works.

Segmentation and Reference Objects
As is well known, the first step in the OBIA approach involves image segmentation to produce homogeneous and discrete objects.Later, these objects, rather than pixels, are used as the classification unit [36].It is important to note that automatic segmentation is a crucial process in OBIA approaches.This process has already been successfully achieved on plastic greenhouses by Tarantino and Figorito [9] through the multi-resolution segmentation algorithm included in the eCognition Developer software (Trimble, Sunnyvale, California, United States).Working on a time series of ASTER images for detecting outdoor crops, Peña-Barragán et al. [24] used the multi-resolution segmentation algorithm implemented in eCognition by equally weighting the green, red and NIR bands of all images with 15-m spatial resolution.In the same way, Vieira et al. [25] carried out their segmentation process by using eCognition's multi-resolution algorithm on a time series of 30-m GSD Landsat images for the detection and mapping of sugarcane in Brazil.In both cases, the agricultural plots had a much larger size than our greenhouses.Therefore, in the present work, the segmentation process will be based on the VHR orthoimage from WV2.
On the other hand, we have to take into account that the final classification accuracy for horticultural crops under greenhouses is going to be highly dependent on the automatic segmentation process carried out on the WV2 orthoimage.Hence, we decided to conduct a manual digitizing in order to obtain the best possible segmentation on our objects of interest.This was already done in a similar way by Aguilar et al. [10].Therefore, 694 reference objects provided by Vicasol were manually digitized over the WV2 pan-sharpened orthoimage with 0.4-m GSD, covering a total of 571 ha (Figure 2).The individual average area of the monitored greenhouses was 8200 m 2 .According to the ground truth data, the existing crops were tomato Solanum lycopersicum L. (54% of the area), pepper Capsicum annuum L. (27%), aubergine Solanum melongena L. (10%) and cucumber Cucumis sativus L. (9%).It is worth noting that tomato and pepper crops were composed of different varieties (e.g., Californian red pepper, Italian green pepper, yellow pepper, orange pepper, plum tomato, cocktail tomato, cherry tomato, truss tomato).

Definition and Extraction of Features
Two different information sources have been tested in this work: Landsat 8 time series and a WV2 single image.Several groups of input variables or features were retrieved from both data sources.These object-based features involved spectral information (mean and standard deviation values of all pixels within an object), several vegetation indices (VIs) related to crop vegetation vigor and textural information based on the gray-level co-occurrence matrix (GLCM) proposed by Haralick et al. [37].In the case of the Landsat 8 time series, the considered object-based features (Table 3) were computed from the seven bands in both visible and SWIR regions included in the eight pan-sharpened Landsat 8 orthoimages.For WV2, the eight bands provided by the MS single orthoimage were used to obtain the object-based features shown in Table 4.In regard to GLCM features, only three out of 14 originally proposed by Haralick et al. [37] were finally considered.Homogeneity, dissimilarity and entropy texture features were always computed considering all of the directions over the NIR region (i.e., NIR1 and NIR2 bands in the case of WV2 and the NIR band when the Landsat time series were analyzed).

Tested Features Description Reference
Indices MSAVI1 (Modified Soil-Adjusted Vegetation Index) NDWI (Normalized Difference Water I) (C−NIR2)/(C + NIR2) [46] Texture GLCMh1 GLCM homogeneity sum of all directions from NIR1 and NIR2 [37] GLCMh2 [37] GLCMd1 GLCM dissimilarity sum of all directions from NIR1 and NIR2 [37] GLCMd2 [37] GLCMe1 GLCM entropy sum of all directions from NIR1 and NIR2 [37] GLCMe2 [37] The OBIA software eCognition v. 8.8 was employed for the extraction of the object-based features from both Landsat 8 and WV2 orthoimages shown in Tables 3 and 4 [37][38][39][40][41][42][43][44][45][46].To this end, the chessboard segmentation algorithm included in eCognition was applied on a thematic layer composed of a previously digitized vector file with the 694 reference greenhouses.Using this approach, the software only projects the vector file on the images in order to obtain a boundary adapted to the pixels that compose them (Figure 3).One eCognition project was conducted for the eight multi-temporal Landsat 8 pan-sharpened orthoimages by halving four times the original 15-m GSD on the created project tab in order to increase the geometric resolution to 0.9375 m per pixel.Thus, the segmentation fit better to the reference greenhouses (see Figure 3a,c,d).In the case of the single WV2 MS orthoimage, another project was performed within eCognition, where the original GSD (1.6 m) was enhanced up to 0.8 m per pixel (Figure 3b).The main goal was always to attain the best possible segmentation for our 694 reference greenhouses.

Decision Tree Modeling and Evaluation
A decision tree (DT) classifier based on the algorithm proposed by Breiman et al. [47] has been used in this work.The major reason for its choice lies in that the results from this classifier are simple to understand.In fact, in Peña-Barragán et al. [24], the DT classifier is described as a "white box", because both its structure (formed by several splits) and its final results (terminal leaves) are easy to interpret.In fact, the DT classifier could be implemented directly in eCognition by means of rule sets.DT is a non-parametric rule-based method in which each node makes a binary decision that separates either one class or some of the classes from the remaining classes.Classification through the DT algorithm is increasingly applied in remote sensing data.In fact, it has been the selected classifier in recent OBIA investigations focused on outdoor crop identification [24,25,48].Other advantages of DT classifiers are: (1) they fit well in the OBIA procedure; (2) DTs are computationally fast and make no assumptions regarding the distribution of the data; (3) they are able to take numerous input variables and perform rapid classification without being severely affected by the "curse of dimensionality"; and (4) the DT methods provide quantitative measurements of each variable's relative contribution to the classification results, so allowing users to rank the importance of input variables.
Once the object-based features corresponding to the 694 reference greenhouses were computed from Landsat 8 time series and WV2 imagery, they were exported by means of eCognition software in order to generate the following three different datasets: (1)  The StatSoft STATISTICA 10 software was used for the classification and regression decision tree analysis (StatSoft Inc., Tulsa, OK, United States).Briefly, the DT algorithm seeks to split the data into segments as homogeneous as possible with respect to the dependent or response variable.In this case, the categorical dependent variable was the kind of crop that had been grown under plastic coverings in the reference greenhouses during autumn 2013 (Table 1).The DT node-splitting rule was the Gini index, which is a measure of impurity for a given node [49].Its application attempts to maximize the homogeneity of the child nodes with respect to the values of the dependent variable.
Cross-validation provided an objective measure of quality for the computed DT model [25,50,51].The experimental design was implemented within the STATISTICA environment through a stratified 10-fold cross-validation procedure, leading to one confusion matrix for each fold.In 10-fold cross-validation, the data are randomly partitioned into 10 equally-sized samples.Of the 10 samples, a single sample is used to validate the model, and the remaining 9 samples are used to train the model.The cross-validation process is then repeated 10 times, with each of the 10 samples used exactly once as the validation data.A final confusion matrix was obtained by summing all of them.The classification accuracy assessment in this work was finally based on this final error matrix provided by cross-validation [52].Hence, the accuracy measures computed were user's accuracy (UA), producer's accuracy (PA) and overall accuracy (OA).Finally, the Fβ measure [53], which provides a way of combining UA and PA into a single measure, was also computed according to Equation (1), where the parameter β determines the weight given to the accuracy computed as PA or UA.The value used in this study (β = 1) weighs UA equal to PA. β 1 β (1)

Classification Accuracy with Regard to Data Sources and Features
Three main datasets were tested in this work (i.e., only Landsat 8 time series data, only WV2 single orthoimage data and, finally, Landsat 8 time series plus WV2 single orthoimage data).In addition, a special case involving only the single Landsat 8 orthoimage taken on 12 September was conducted just to be compared with the single WV2 case taken also in September.Moreover, four strategies or groups of features were considered for each dataset: (1) all features (All); (2) all features without GLCMs (no GLCM); (3) all features without GLCMs, mean values (no GLCM, mean); and (4) all features without GLCMs, mean and standard deviation values (no GLCM, mean, SD).
Figure 4 shows the accuracy per class and OA based on the 10-fold cross-validation method for the different datasets tested.It can be noted that the best OA values were attained with the two first strategies.Furthermore, the second-order textural parameters based on GLCM added no information in order to improve the classification accuracy.Working on outdoor crops identification, Peña-Barragán et al. [51] reported that texture features achieved extremely poor overall accuracy as compared to spectral features.However other research found that texture attributes based on GLCM were relevant to discriminate between some crops with similar spectral responses [24,25].In this sense, it is important to bear in mind that the calculation of GLCM is computationally intensive and very time consuming.Notice that, apart from GLCM descriptors, textural information is somehow included by means of SD values, which can be considered as first-order textural features.The best OA value (81.3%) was reached when using no GLCM strategy on Landsat 8 time series (Figure 4a).When the Landsat 8 + WV2 dataset was employed, a slightly lower OA value of 79.4% was achieved with the no GLCM strategy (Figure 4b).OA values were far worse when only a single growing season period was considered, highlighting the suitability of multi-temporal images taken along the growing season to improve under greenhouse crop discrimination.In this way, and applying the no GLCM strategy, OA values of 75.1% and 67% were achieved for WV2 and Landsat 8 September datasets, respectively (Figure 4c,d).The better classification results found in the case of WV2 can be attributed to its much higher spatial resolution.Regarding the results from the third strategy (no GLCM, Mean) depicted in Figure 4, it should be noted that object spectral information based on mean values could not be removed without causing a decrease in OA close to 2%.This result strengthens the importance of applying proper atmospheric corrections on the original images.The OA values of around 80% achieved in this work are rather promising as compared to those reported by other authors working with a similar methodology, but on outdoor crops.For instance, Peña-Barragán et al. [24] reached an OA of 79% applying the DT classifier for detecting 13 major crops cultivated in the agricultural area of Yolo County, California, USA.This classification accuracy rose up to 88% by using support vector machine (SVM) or multilayer perceptron (MLP) neural network classifiers [51].In another study, twelve outdoor crops and two non-vegetative land use classes were classified with 80.7% OA at cadastral parcel levels through multi-temporal remote sensing images [48].Furthermore, working on mapping outdoor crops (classes corn, soybean and others), an OA of 88% was reported by Zhong et al. [2].Finally, Vieira et al. [25] applied the DT classifier to Landsat time series in a binary classification approach between sugarcane and other crops, assessing an OA of 94%.
Regarding the accuracy per class, Figure 4 shows that tomato and pepper crops were fairly well identified for all datasets.Again, and always using the no GLCM strategy, the best Fβ values for tomato (87.7%) and pepper (87%) were attained for the Landsat 8 time series dataset (Figure 4a).With the same strategy, the Fβ measure reached values of 85.1% and 84.5% for tomato and pepper, respectively, when using the Landsat 8 time series plus WV2 dataset.It is important to underline that Fβ values for tomato and pepper crops were significantly worse for the datasets where only a single image of WV2 (80.3% tomato and 81.5% pepper) or Landsat 8 September (76.2% tomato and 68.8% pepper) were employed.On the other hand, the other two crops presented much poorer classification accuracies.The best Fβ value for aubergine identification (60.3%) was achieved using the WV2 dataset with no GLCM strategy.Moreover, an accuracy of 63.2% in terms of Fβ was reached in the case of cucumber for the Landsat 8 time series dataset and, again, through the no GLCM strategy.Working on 13 outdoor crops, Peña-Barragán et al. [24] reported Fβ values ranging from 94% for rice to 66.7% for rye, achieving an accuracy of 80% in the case of tomato.
More detailed information about the accuracy per class is presented in Table 5, showing the corresponding confusion matrix for the Landsat 8 time series dataset and no GLCM strategy.It can be noted that a large number of under-greenhouse aubergine crops are misclassified as tomato.In fact, aubergine and tomato crops have a similar transplant date in Almería, usually in August (Tables 1 and 2).Both crops are quite heat tolerant, and in the study area, greenhouse coverings only need a slight white paint layer to protect the crops during the transplant stage.The white paint is usually washed in September or October.Regarding cucumber crops, they have very variable transplant dates, ranging from August to October (Table 2).For early transplanted cucumber crops, the white paint is more intense to protect the small plants during their first stage.A very slight white paint is usually applied when the cucumber crops are planted in early to mid-September and even not necessary when the cucumber plants are planted very late.Furthermore, in the study area and regarding cucumber crops, it is common to rely on a double-roof in the greenhouse, which reduces radiation losses.All of these actions cause a highly variable spectral signature on the shape of the time series profile recorded from greenhouse coverings.Finally, pepper crops are mainly planted during July and August.This is a crop that needs a significant quantity of white paint to enhance the growth of the plant, as well as to protect the transplant stage.Without a doubt, this is the most unique characteristic for the autumn pepper crop in the "Poniente" region of Almería.The pepper crops are mainly misclassified as tomato or cucumber (Table 5).In Figure 5 is depicted the evolution of the time series profile.The cucumber greenhouse located in the center, planted on 27 August 2013, appears without white paint on 11 August (Figure 5b,f).It seems to be whitewashed in Figures 5a,c,e,g (30 and 12 September, respectively), and finally, this greenhouse covering is shown without any trace of paint in November (Figure 5d,h).Moreover, November (Figure 5h) turns out to be the month when most of the greenhouses allow one to visually perceive (false color images) the presence of crops through an increase in the NIR reflectance.This is not the case of the greenhouse located just on the right of the aforementioned cucumber one.In fact, it is a tomato crop that was whitewashed in September, and it had not yet been washed in November (Figure 5h).
The capacity of the DT approach was used to select the object-based features and their cutting thresholds that best fit to every type of under-greenhouse crop.In view of the complexity of the computed DT, presenting up to 34 terminal nodes for the case shown in Table 5 (Landsat 8 time series dataset and no GLCM strategy), just the first divisions of that DT (only nine terminal nodes) are presented (Figure 6).At the bottom and outside of each black box in Figure 6 appears the object-based feature that is causing the next split in the decision tree.Below, the readers can see the figure that represents the cutting threshold for the split.Oval nodes are terminal nodes belonging to a particular class.Thereby, the OA value attained with the pruned DT shown in Figure 6 was 75.9%, considering the last box (N = 419) as the Class 4 (Tomato) terminal node.It is worth noting that the object-based feature located at the top of the tree was the mean value of SWIR2 for the Landsat 8 pan-sharpened orthoimage taken on 11 August (AGO11 Mean SWIR2).SWIR2 reflectance values higher than 44.85% (right branch of the tree) were mainly associated with greenhouses that received an important dose of white paint in August, so pointing to pepper crops.In fact, 149 out of the 191 selected greenhouses with higher SWIR2 turned out to be pepper crops.The remaining high SWIR2 objects in August were identified as tomato (33 greenhouses belonging to Class 4), aubergine (eight greenhouses for Class 1) and cucumber (one greenhouse belonging to Class 2).The next steps in that branch of the DT were addressed to separate these tomato greenhouses.In this way, features, such as the mean value for the green band reflectance on 24 June, the Green Vegetation Index (GVI) on 26 July and GVI in November, were chosen to carry out this task.At the end of June, pepper greenhouses were being prepared to receive the plants, so they were usually bare and presented a lower value for the green band reflectance.Regarding GVI values for the pepper crops, they were often higher at the end of July and lower in November as compared to other crops in the study area.On the left branch of the DT (Figure 6) we observe most of the aubergine (72), cucumber (61) and tomato (327) crops, which received low doses of white paint (SWIR2 reflectance ≤ 44.85%), but it also contained 43 greenhouses of Class 3 (Pepper).The SWIR1 mean value on 26 July was used to extract these pepper crops with a high dose of white paint when transplanted early.The rest of the object-based features used in this branch (i.e., GNDVI in November, mean NIR in September and Normalized Differential Senescent Vegetation Index (NDSVI) in May) were aimed at discriminating tomato, cucumber and aubergine crops.The last box in Figure 6 still contains 419 objects classified as tomato in the DT with nine terminal nodes.The whole DT, presenting 34 terminal nodes, continues from that point using the other 23 features extracted along different periods, such as mean reflectance values (NIR, Blue, SWIR2 and SWIR1), standard deviations (SWIR2, blue, green and NIR) and VIs (NDSVI, Normalized Difference Tillage Index (NDTI), GVI and GNDVI).
In view of the above, the identification of under-greenhouse pepper crops can be carried out with high success rates, mainly based on the crop management techniques applied during the growing cycle for each agricultural region.In our test field, an OA of 92.7% and a kappa coefficient of 0.81 were achieved considering only two classes (pepper and others) by using Landsat 8 time series with no GLCM strategy.In addition, a very simple DT (similar to the top part of Figure 6) with only six terminal nodes and only five features (i.e., AGO11 mean SWIR2, JUL26 mean SWIR1, JUN24 mean GREEN, SEP12 SD BLUE and JUL26 GVI) were required for the achievement of the aforementioned results.

Optimal Time Period for Collecting Satellite Images
This section is devoted to detecting the optimal temporal windows for gathering the satellite imagery time series, as well as crucial time periods for under-greenhouse crop identification.The first step was to assess the classification accuracy attained from each of the nine available single orthoimages, always using the no GLCM strategy.Table 6 shows that the best classification results were obtained from the WV2 orthoimage, likely due to its much higher spatial resolution.With regard to the eight Landsat 8 orthoimages, the best OA values were assessed from the images taken in August (73.2%),on 26 July (69.7%) and in September (67%).The OA values for the rest of the images ranged from 65.4% to 63.2%.According to the agricultural practices of the working area, the importance of the images taken in August and late July could lie in the application of white paint.Along September, some crops reach a significant size, and therefore, leaf pigments (chlorophyll, carotenoids and anthocyanins) begin to influence the signal captured by the satellites (e.g., NIR band), although in this month, there are still numerous greenhouse coverings painted white (Figure 5g).In November, most of the greenhouses have been washed, so being transparent to the electromagnetic radiation reflected from the crops inside (Figure 5h).However the OA value computed from this single period was quite poor.Table 7 shows the accuracy results for different multi-temporal satellite imagery combinations.The best combination of three Landsat images included the scenes corresponding to 26 July, August and September, achieving an OA of 75.5%.The time series profile was improved by adding the November information, the OA value increasing to 77.8%.The information contained on the scene taken in June, the period related to the greenhouse preparation for the transplant stage, also turned out to be relevant to improve OA up to a value of 79.8%.Information about previous crops given by the May scene did not raise the final accuracy.When the WV2 dataset was also added to Landsat time series, the best OA values (81.1% and 82.3%) were achieved.Likely, including more temporal data within the growing season from June 2013 to March 2014 would have led to improving the current results.It is worth noting that a couple of medium spatial resolution super-spectral satellites, named Sentinel-2, is going to be launched in 2015-2016 by the European Space Agency.These satellites will be able to acquire images of the Earth covering 13 spectral bands with spatial resolutions of 10 m (four visible and NIR bands), 20 m (six red-edge/SWIR bands) and 60 m (three atmospheric correction bands).The revisit time will be only five days from the two-satellite constellation with open and free data access [54].
Another issue to consider in further works would be the use of seasonal statistics [55] or some automated crop mapping algorithms based on the shape of the time series profile of VIs [2], which should be more consistent over years than snapshot values the original images.Both options could contribute to making the multi-temporal approach less dependent on the availability of particular acquisition dates and less susceptible to missing data on crucial dates (e.g., due to cloud presence or incomplete spatial coverage).

Conclusions
As far as we know, this work is the first attempt to identify horticultural crops that are growing under greenhouse plastic coverings by using multi-temporal and multi-source satellite images.To this end, several spectral, textural and VIs object-based features were extracted from a WV2 orthoimage and eight multi-temporal Landsat 8 pan-sharpened orthoimages by using an OBIA approach.These previously-segmented objects were subsequently analyzed by means of the DT classifier.When the full Landsat 8 time series was applied, classification results reached an OA of 81.3% from performing the 10-fold cross-validation method over 694 reference greenhouses or the ground truth.The classification results were quite good in the case of tomato and pepper crops, presenting Fβ values of around 87%, but worse on cucumber and aubergine crops, with Fβ values ranging from 50% to 63%.These results hardly improved by adding the information of the WV2 single image, which had been used to undertake the initial segmentation stage.In fact, the best OA reached 82.3% (kappa coefficient of 0.71) when the WV2 data were added to five Landsat images (8 June, 26 July, August, September and November).The most relevant object-based features were related to spectral information (mean reflectance and SD values) and VIs.The most important spectral information for correct classification of different crops under greenhouses was mainly due to the greenhouse management practices.As regard VIs features, GVI, GNDVI and NDSVI turned out to be the most significant ones.The second-order textural parameters based on GLCM added no information in order to improve the classification accuracy.
As is widely known, outdoor crop classification is rather complex due to agronomic factors and farming decisions.Different crops can exhibit very similar developmental patterns and growth calendars.In addition, the same crop may be sown on different dates due to the farmers' decision.All of these variables become much more complex in the case of crops under greenhouses, where the existence of plastic coverings and their particular management during the growing season make the crop identification via remote sensing more difficult.The most defining characteristic for crop identification in this work turned out to be the unique whitewashing management in the case of pepper crop.Thus, high reflectance values in the SWIR1 and SWIR2 bands in the months of July and August were related to the white paint applied to the plastic coverings of this crop.In this regard, pepper could be discriminated from other crops through a DT binary classification by using Landsat 8 time series and only five features (AGO11 mean SWIR2, JUL26 mean SWIR1, JUN24 mean green, SEP12 SD blue and JUL26 GVI), attaining an OA of 92.7% and a kappa coefficient of 0.81.
Although it has been proven that satellite image resolution was relevant to the final results (i.e., the best OA on single images was achieved using the WV2 dataset), the use of multi-temporal images along growing season turned out to be crucial to improve the discrimination between different crops.In that sense, the use of pan-sharpened orthoimages obtained from free available Landsat 8 scenes resulted in being sufficiently successful.However, bearing in mind the small area of plastic greenhouses and the complexity of the task addressed in this work, it seems obvious that using higher resolution time series, both temporal and spatial, would result in better outcomes.In that way, the new Sentinel-2 satellites could offer a valuable source of data in the near future.
The optimal temporal windows for collecting the satellite imagery time series in order to identify the autumn greenhouse crops should include more temporal data within the growing season from June to March.It is important to note that greenhouse coverings are usually whitewashed again in February.Regarding our imagery time series from May to November, the crucial collection periods proved to be July and August.Likely, adding more satellite images within this period of time could lead to better results.
The methodology applied in this work could be adapted to other intensive agricultural areas using plastic covering, such as Morocco (e.g., Bir Jdid, Oulidia, Essaouira and Agadir) or China.However, further works are needed to better understand the influence of the type of plastic (material, color, thickness) on crop detection.

Figure 1 .
Figure 1.Location of the study area on a RGB Landsat 8 image taken in August 2013, covering the "Poniente" region.Coordinate system: ETRS89 UTM Zone 30N.

Figure 4 .
Figure 4. Per-class accuracy given by the corresponding Fβ index and overall accuracy (OA) using different groups of features and data sources.(a) Landsat 8 time series (eight orthoimages); (b) Landsat 8 time series plus WV2 single orthoimage; (c) WV2 single orthoimage taken on September 30; (d) Image Landsat 8 taken on September 12.

Table 1 .
Characteristics of the 694 reference greenhouses.SD: Standard deviation.

Table 2 .
Schematic summary of the typical transplant pattern for the four studied crops.

Table 5 .
Confusion matrix for Landsat 8 time series dataset and no GLCM strategy.PA, producer's accuracy.

Table 6 .
Accuracy indices from single orthoimages using the no GLCM strategy.

Table 7 .
Accuracy indices from different combinations of multi-temporal orthoimages using the no GLCM strategy.