1. Introduction
The extensive and steadily expanding use of plastic films in agriculture, and particularly in protected horticulture, is reported worldwide since the middle of the twentieth century [1]. Due to their synoptic acquisitions and high revisit frequency, the data obtained by remote sensing can offer a significant contribution to provide periodic and accurate pictures of the agricultural sector [2]. Therefore, an accurate remote sensing mapping of plastic covered greenhouses (PCG), walk-in tunnels, low tunnels and plastic-mulched farmland (PMF) has important practical significance for analyzing land-use evolutions and guiding regional agricultural production [3,4].
The recent work published by [5] provides an overview of the general research dynamics regarding the topic of remote sensing-based mapping of agricultural greenhouses and plastic-mulched crops throughout the 21st century. Several research lines have emerged throughout the 21st century. They are based on (i) input satellite data sources, including optical (e.g., Landsat or Sentinel-2) and SAR (Radarsat-2, Sentinel-1)imagery; (ii) image processing approaches (i.e., pixel-based or object-based image analysis (OBIA)); (iii) input vector features for image classification (multispectral or hyperspectral, geometrical, texture, 3D information, radar, lidar and data fusion); (iv) image classification methods (unsupervised or supervised, parametric or non-parametric and machine learning or deep learning-based approaches); (v) single or multitemporal remotely sensed data; vi) and only mapping plastic covering structures or going a step further to classify and map plastic covered crops.
A groundbreaking very high resolution (VHR) satellite named WorldView-3 (WV3) was successfully launched on 13 August 2014. WV3 opened up new opportunities for analyzing and mapping plastic materials thanks to its global super-spectral observation capabilities (eight VNIR bands and eight SWIR bands) [6].
Throughout the last decade, studies on plastic detection via remote sensing have been focused on anthropogenic marine debris mapping in the oceans [7,8,9]. Acuña-Ruz et al. [7] employed eight VNIR bands and the panchromatic band of WV3 data, along with different classification methods, to map plastic debris on beaches. Biermann et al. [8] used remote sensing techniques with Sentinel-2 for mapping patches of floating macroplastics on the ocean surface. Themistocleous et al. [9] detected plastic targets on the sea surface through the proposed Plastic Index (PI) based on Sentinel-2 data. Moreover, several works have been carried out for mapping plastic films in agriculture, such as PCG [10,11] or PMF [12,13]. Novelli et al. [10] carried out the first work based on Sentinel-2 Multispectral Instrument (MSI) and Landsat 8 OLI images to map PCG by adopting an OBIA approach and an RF classifier. Yang et al. [11] developed a specific index for the detection of PCG (Plastic Greenhouse Index; PGI) from Landsat TM and Landsat 8 OLI images. Lu et al. [12] developed a simple but robust decision tree (DT) classifier to detect PMF from only two Landsat 5 TM images. Hasituya et al. [13] used Landsat 8 OLI imagery and GoogleEarth data (training phase) to map PMF using spectral and textural characteristics. Research by [14] addressed the classification of plastic materials in urban areas using the short-wave infrared (SWIR) bands from a WV3 image, obtaining promising results. Asadzadeh and de Souza Filho [15] tested the WV3’s SWIR bands for the detection of hydrocarbons, which are spectrally similar to plastics. In fact, plastic is a type of synthetic polymer or hydrocarbon commonly made by linking numerous hydrocarbon monomers together into long chains of molecules [16].
Several indices have been proposed in the aforementioned studies for mapping plastic materials. Among them, we highlight the following: (i) The Plastic Greenhouse Index (PGI), derived from the visible and near-infrared (VNIR) bands of Landsat 7 Enhanced Thematic Mapper plus (ETM+) imagery [11]; (ii) the Plastic Index (PI), from the VNIR bands of Sentinel-2 [9]; (iii) the Plastic–Mulched Landcover Index (PMLI), from the SWIR bands of Landsat 5 Thematic Mapper (TM) images [12]; (iv) the Relative-absorption Band Depth index (RBD), derived from WV3 SWIR bands [15]; and (v) the Normalized Difference Plastic Index (NDPI), from WV3 SWIR bands [14].
Regarding the spectral properties (350–2500 nm) of human-made plastic materials used in covered agriculture, Levin et al. [17] analyzed, in the laboratory, 15 samples of polyethylene (PE) sheets (transparent, black, opaque, and yellow opaque), and various nets (white and black) used in Israel, utilizing a field spectroradiometer Analytical Spectral Devices FieldSpec Pro (ASD 2001). They also considered field data taken from a hyperspectral AISA-ES image with a spatial resolution of 1 m. They pointed out three major absorption features around 1218 nm, 1732 nm and 2313 nm for plastic detection. Moreover, these valley values were not affected by settling dust, whitewashing or by the underlying surface. They marked the value of 1732 nm as the best spectral feature for plastic mapping.
More recently, Guo and Li [14] analyzed the spectra corresponding to urban plastic materials. They found that the reflectance of different plastics was highly variable in the visible spectral range depending on, for instance, the colors of the plastic samples. However, in the NIR to SWIR spectral range, the shape of spectral curve was related to the specific plastic composition. Plastics are organic compounds with numerous C-H bonds. Existing studies have indicated that these bonds have diagnostic absorption features at the wavelength regions of 1100–1250 nm (second overtone of C-H stretching mode), 1300–1450 nm (C-H combination band), 1600–1800 nm (first overtone of C-H stretching mode), and 2150–2500 nm (combination band) [18]. Guo and Li [14] tried to transfer this knowledge to the spectral information included in the eight WV3 SWIR bands, detecting two declining trends: (i) between band 10 (SWIR10) and band 12 (SWIR12) (i.e., 1570 nm to 1730 nm), and (ii) from band 13 (SWIR13) to band 16 (SWIR16) (i.e., 2165 nm to 2330 nm). These findings led to the proposal of the NDPI index [14].
The innovative goal addressed in this work relies on testing the predictive power attributed to features derived from the WV3 SWIR bands, in contrast to those obtained from the WV3 VNIR bands for mapping PCG within the context of an OBIA approach. Thus far, and to the best knowledge of the authors, this is the first research work that deals with this topic.
2. Study Area
The study area, located in the province of Almeria (Southern Spain), comprises a rectangle of 2349.6 m per 8787.6 m, covering an area of 20.65 km2 centered on the geographic coordinates (WGS84) 36.7774°N and 2.6572°W (Figure 1). This study area was previously selected in different research projects related to PCG mapping because it is just at the core of a great concentration of greenhouses in Almeria. In fact, more than 70% of the area is covered by PCG, though there are also roads, small cities, bare soil, irrigation ponds and vegetation.
3. Datasets and Pre-Processing
3.1. WorldView-3 Image
A single WorldView-3 (WV3) bundle image (PAN + VNIR + SWIR) in Ortho Ready Standard Level-2A (ORS2A) format was acquired on June 11, 2020. This sensor has one panchromatic band (PAN) with a resolution of 0.31 m (at nadir), eight visible and near infrared bands (VNIR) with a resolution of 1.24 m and eight shortwave infrared bands (SWIR) with a resolution of 3.7 m (Table 1). The off-nadir angle of the WV3 image was 27.6°, with collected pixel sizes of 0.381 m (PAN), 1.522 m (VNIR) and 4.585 m (SWIR), while the final product presented pixel sizes of 0.3 m, 1.2 m and 3.7 m for PAN, VNIR and SWIR, respectively.
From the PAN and VNIR WV3 images, a pansharpened image with 0.3 m ground sample distance (GSD) was attained by means of the PANSHARP module included in Geomatica v. 2018 (PCI Geomatics, Richmond Hill, ON, Canada). For the three WV3 images (i.e., pansharpened, VNIR and SWIR), the coordinates of seven ground control points (GCPs), obtained by differential global positioning systems, were used 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 Andalusian Government, was used to carry out the orthorectification process. Three WV3 orthoimages with 0.3 m GSD (pansharpened), 1.2 m (VNIR) and 3.7 m (SWIR) were generated using Geomatica v. 2018. Finally, the VNIR and SWIR orthoimages were atmospherically corrected by using the ATCOR (atmospheric correction) module included in Geomatica v. 2018. This absolute atmospheric correction algorithm involves the conversion of the original raw digital numbers to ground reflectance values by applying the MODTRAN (MODerate resolution atmospheric TRANsmission) radiative transfer code [19].
3.2. Ground Truth
A ground truth or reference dataset was manually digitized from the 0.3 m WV3 pansharpened orthoimage (Figure 1c). It contained 1979 polygons representing individual agricultural greenhouses (GH class), and one polygon, including the non-greenhouse class (Non-GH).
4. Methodology
The methodology proposed in this article addresses the following aspects: (i) Pre-processing protocol for WV3 images (explained in Section 3.1); (ii) segmentation of the WV3 VNIR orthoimage and the ground truth for delineating greenhouses; (iii) binary classification of the entire study area into two classes (GH and Non-GH) by applying an OBIA approach based on the VNIR and SWIR WV3 orthoimages. A flowchart describing the methodology proposed is shown in Figure 2.
4.1. Segmentation
The first stage in an OBIA approach involves image segmentation to produce homogeneous objects which will be used afterwards as classification units [20]. In previous works dealing with mapping PCG, this segmentation was performed on atmospherically corrected VNIR VHR satellite imagery [21,22] using the widely known multi-resolution (MRS) image segmentation algorithm implemented in Trimble eCognition Developer v. 9.5 (Trimble, Sunnyvale, CA, USA). Bearing in mind, that the final classification accuracy for PCG is going to be highly dependent on the segmentation goodness. We decided to use the manually digitized ground truth (Section 3.2) in order to obtain the best possible segmentation for our objects of interest (i.e., PCG). In this way, the 1979 GH polygons that made up the ground truth were uploaded as a vector thematic layer into a project developed into eCognition. After that, the polygons were transferred to the WV3 ATCOR orthoimages (VNIR+SWIR) using the chessboard segmentation algorithm and the thematic layer. In addition, the area belonging to the Non-GH class had to be segmented as well. The adopted solution was to choose the best MRS segmentation attained for the whole study area for the Non-GH area. The outcome of the MRS algorithm, which is a bottom-up region-merging technique starting with one-pixel objects or seeds, is controlled by three main factors: (i) scale parameter (scale), determining the maximum allowed heterogeneity for the resulting segments; (ii) weight of color and shape criteria in the segmentation process (shape); and (iii) weight of the compactness and smoothness criteria (compactness). Thus, these factors and the bands’ combination with their corresponding weights had to be set for this work. According to [23], the best bands combination was fixed at equally weighted blue-green-NIR2 bands from the VNIR ATCOR WV3 orthoimage. Moreover, compactness was fixed at 0.5 [24]. Thousands of segmentations from applying the MRS algorithm were computed by means of a semi-automatic eCognition rule set characterized by a looping process which varied shape (values ranged from 0.1 to 0.5, with a step of 0.1) and scale (from 40 to 600, with a step of 1) tuning parameters. The selection of the best scale and shape parameters for the WV3 orthoimage was carried out through a modified version of Euclidean Distance 2 (ED2), a supervised segmentation quality metric included in a command line tool named AssesSeg [10]. It is worth noting that the ED2 measure, originally proposed by [25], has provided very good results for working on PCG segmentation [21,22]. In this work, and according to [24], we used 400 PCG objects from the manually digitized ground truth as the reference, set to run AssesSeg. Full details about the modified ED2 measure and the standalone command line tool (AssesSeg.exe) can be found in the work published by [10].
Summing up, our idea was to obtain the best MRS segmentation representing greenhouses in the study area, improving it in the case of the GH class by using manually digitized ground truth. This segmentation had 3188 objects.
4.2. Definition and Extraction of Object-Based Features
The object-based features used in this work (Table 2) to cope with the binary OBIA classification in GH and Non-GH classes were retrieved from both the WV3 VNIR and the WV3 SWIR orthoimages through the aforementioned eCognition project. Note that, in that project, the GSD for all 16 bands was automatically reset to 1.2 m GSD (the smallest GSD for VNIR and SWIR orthoimages).
The basic spectral information (33 object features) consisted of the mean and standard deviation (SD) values of the ground reflectance values (ranging from 1 to 10,000) for all the pixels within an object for each of the 16 bands, and the VNIR brightness (average reflectance values for all the eight VNIR bands) was extracted from the eCognition project. Additionally, six object features consisting of plastic and vegetation indices were also derived from the WV3 orthoimages. Table 2 shows all 39 object features tested in this work, differentiating whether they were extracted from the VNIR or SWIR orthoimages.
4.3. Decision Tree Modeling and Classification Accuracy Assessment
A supervised decision tree (DT) classifier based on the algorithm proposed by [28] was used for the three tested strategies: (i) using only the 20 object features derived from VNIR; (ii) using only the 19 object features derived from SWIR data; and (iii) using all 39 features (VNIR+SWIR) (Table 2). Although other supervised classifiers such as random forest or support vector machine (black box classifiers) could achieve better classification results, DT provides clear decision rules with fixed threshold values that can be used in further works without any training phase. The supervised classification through the DT machine-learning method has already been applied in PCG [21] and PMF mapping [12], with good and robust results.
STATISTICA v10 (StatSoft Inc., Tulsa, OK, USA) was the software used for the classification decision tree analysis. Aguilar el al. [21] and Nemmaoui et al. [22] already used this tool for dealing with PCG detection. The DT node-splitting rule was the Gini index, which is a measure of impurity for a given node [29].
From the 3188 objects attained in Section 4.1, a balanced sample containing 1175 segments for the GH class and 1175 for Non-GH class was selected. After that, these objects were classified according to the information of the ground truth (Figure 1c). All these objects represented meaningful objects for each class avoiding mixed segments. The selected Non-GH objects corresponded to bare soil, buildings, vegetation and roads.
The selected dataset (1175 + 1175 objects), with its extracted object features, was used to build the optimal DT for each strategy into STATISTICA v10. In this way, the categorical dependent variable for the DT model was the class (GH or Non-GH), while the object features were selected as the continuous predictor variables. A stratified 10-fold cross-validation procedure, where all of the data are randomly partitioned into 10 equally sized samples, was carried out [3,30]. Any object was used in this stage for validation. The DT classifier provides both a confusion matrix and an assessment of the relative importance of the different features in the classification process.
Since this confusion matrix was based on meaningful objects, the errors related to the segmentation process and/or mixed objects were not measured. Moreover, all 2350 objects were considered for the training phase. That is why, according to [21,22], a new and more realistic classification accuracy assessment based on the real ground truth comprising the entire study area (Figure 1c) was carried out. In this accuracy assessment, the DT rules were implemented in eCognition to classify the 3188 objects in GH or Non-GH classes, following an object-based classification. Then, the vector ground truth was finally exported as a raster file with 1.2 m pixel size. This raster file was used into eCognition to attain a pixel-based error matrix based on the TTA (test and training area) Mask. Overall accuracy (OA) and kappa coefficient (kappa) [31] were computed in both the object-based and pixel-based accuracy assessment tests. Additionally, user’s accuracy (UA) and producer’s accuracy (PA) for the GH class were also computed for the pixel-based accuracy assessment [31].
5. Results
5.1. Segmentation
The scale, shape and compactness parameters set to perform the MRS segmentation were 403, 0.5 and 0.5, respectively, yielding a modified ED2 value of 0.313 and a segmentation composed of 3188 objects for the entire study area. Figure 3a depicts a detailed view of this MRS segmentation. Considering that we had the ground truth (Figure 3b), we decided to improve the original MRS segmentation of GH objects through vector geometry coming from ground truth to obtain the final segmentation (Figure 3c), where the GH polygons fully matched the digitized ones in the ground truth. Once the GH objects were segmented, the Non-GH area (green area in Figure 3b) was segmented by applying the MRS with the aforementioned scale, shape and compactness parameters. Note that since the area to be segmented, in this case, was not exactly the same as in the initial segmentation (i.e., segmentation-only covered Non-GH area), the results also varied slightly.
5.2. Object-Based Accuracy Assessment
Table 3 shows the confusion matrices from the DT classifier, computed through the 2350 meaningful objects (1175 GH and 1175 Non-GH). It is clear that the best strategy was to use all the information included in the VNIR and SWIR bands of WV3 (i.e., All Features), obtaining values of 96.64% and 0.933 for OA and kappa, respectively. Slightly worse accuracies were achieved (OA = 95.06% and kappa = 0.901) when only SWIR features (i.e., SWIR strategy) were applied. In the case of the VNIR strategy, the accuracy values plummeted to 89.62% and 0.792 for OA and kappa, respectively.
5.3. Importance of Features and Decision Tree Models
The rankings of the 10 most important features used for the three tested strategies are depicted in Table 4. NDPI was the most relevant index for both All Features and SWIR strategies, closely followed by PMLI. However, there was not a third important index in the SWIR strategy, while the All Features strategy had Mean C, PGI, Mean B and even Brightness, as well-ranked features. In the case of VNIR strategy, the SWIR derived features, such as NDPI and PMLI, were replaced by Mean C, Mean B and PI, showing the importance of coastal blue and blue visible bands for PCG mapping.
The DT models developed for each of the strategies had five (SWIR) or six (VNIR and All Features) terminal nodes. NDPI, PLMI and PGI were the features located in the first splits (i.e., the most important) of the DT in the All Features strategy (Figure 4). NDPI and PLMI materialized the first splits for the SWIR strategy, while Mean C and Mean B were in the main places for the VNIR strategy.
5.4. Pixel-Based Accuracy Assessment
A final accuracy assessment was carried out over the whole study area using the DT models attained in the previous section to achieve the object-based classification. At this time, accuracy measures were computed on pixels rather than objects.
Table 5 shows the confusion matrices for this accuracy assessment. The results were slightly better than those shown in Table 3. It is important to note that the OA and kappa values in Table 5 are not only statistical estimators but real values, since they are calculated from the entire study area and not from a sample of it.
Again, the best strategies were those that included information from the WV3 SWIR bands (i.e., All Features and SWIR), with the NDPI index being the first and most important cutoff value to classify GH and Non-GH. Moreover, the UA and PA values for the GH class (GH UA, GH PA) resulted in fairly balanced values.
Figure 5 shows the capabilities of the NDPI index to classify plastic materials. All the PCG in Figure 5a presented an NDPI value ranging from 0.225 to −0.035 (Figure 5b). It is important to note that, in June, there were already some whitewashed greenhouses receiving pepper and tomato plants during the months of July and August. In these whitewashed greenhouses, the spectral signature of the plastic cover is masked and, therefore, the value of the NDPI drops significantly. In this way, the whitewashed greenhouses in Figure 5a presented a lower NDPI value (shades farther from the white color) in Figure 5b. The two farm buildings (non-plastic material) marked with the blue letter B in Figure 5a had NDVI values lower than -0.11. In addition to greenhouses, other objects presented high NDPI values in Figure 5b. Among them, we can see: (i) an agricultural building with a plastic cover (Figure 5c); (ii) a local plastics recycling plant (Figure 5d); (iii) piles of agricultural plastic crates (Figure 5e, white ellipse); (iv) plastic solar panels where plastic is used for connecting components, including thrust washers, electrical insulators, pipes, valves and other fittings (Figure 5e, yellow ellipse); or (v) an abandoned greenhouse with its plastic cover very run down and broken (Figure 5f).
After analyzing the classification results, it seems clear that the NDPI index based on the WV3 SWIR bands worked perfectly in the detection of plastic materials. However, many misclassification errors were attained by applying the DT model developed for the best strategy (All Features) depicted in Figure 4. Figure 6a shows the classification results of the All Features strategy in a detailed area of 2000 m per 2000 m, where seven GH were misclassified as Non-GH (yellow color). These errors were mainly due to the masking effect because of the agricultural practice of greenhouse whitewash (greenhouse shading) in some greenhouses.
Although it is not the objective of this work, we decided to propose a new cutoff rule based on an index (NDPI_B) composed of NDPI and Brightness (Equation (1)) to solve the whitewash masking effect.
NDPI_B = NDPI + (Brightness/30000)(1)
In those whitewashed greenhouses, where the plastic covering material is masked (low NDPI values) but presents an important brightness (high Brightness values), the NDPI_B index considers the values of the NDPI and Brightness to improve the classification. Objects with NDPI_B values greater than 0.036 are classified as GH, with the rest belonging to the Non-GH class. Figure 6b depicts the classification results following this knowledge-based DT. As compared with Figure 6a, there were a few more misclassification errors in purple color (Non-GH objects classified as GH) and less in yellow color (GH objects classified as Non-GH). Between the purple group, we have a few of the cases shown in Figure 5 and three big ponds covered by plastic. In yellow, only one greenhouse is shown; after visiting it, it was found to be made of glass and not plastic.
Table 6 shows the confusion matrix for the NDPI_B knowledge-based DT. Comparing with All Features strategy, the OA and kappa were slightly improved. However, the main improvement was related to the producer’s accuracy or error of omission for GH class (GH PA), which represents how well reference pixels of the GH class from the ground truth are classified. This accuracy measure was 96.90% for the All Features strategy (Table 5), rising up to 99.28% using NDPI_B (Table 6).
Finally, Figure 7 shows the reflectance values (mean object figures) for the eight SWIR bands of WV3 image for four PCGs (objects). The two declining trends reported by [14] between SWIR10 and SWIR12, and from SWIR13 to SWIR16, for the three plastic films with different colors (grey, green and black), but without whitewashing, are apparent. These downward trends are not present in the whitewashed PCG (color red in Figure 7). This masking effect produced by whitewash (lime application) is the reason why the NDPI index does not work properly in relation to whitewashed greenhouses.
6. Discussion
The automatic MRS segmentation carried out in this work achieved a modified ED2 value of 0.313 with a resonable visual quality (Figure 3a). In others works, when dealing with VNIR WorldView-2 (WV2) (2 m GSD) or WV3 (1.2 m GSD) atmospherically corrected orthoimages, ED2 values ranged from 0.11 to 0.29. In fact, Aguilar et al. [21], using all eight equally weighted bands and the original ED2 metric proposed by [25], reported the best original ED2 value (0.11) for a VNIR WV2 orthoimage taken in September 2013, and the worst one (0.29) from a VNIR WV2 orthoimage collected in July 2015. Novelli et al. [10] improved the last ED2 value to 0.199, working with the same WV2 orthoimage (July 2015) but using a blue-green-NIR2 equally weighted band combination and the command line tool AssesSeg. Finally, Aguilar et al. [24] reported a modified ED2 valued of 0.141 from a VNIR WV3 orthoimage taken on July 5, 2016; in this case, using the band combination recommended by [10] and AssesSeg. All these data related to the segmentation of PCG reflect that this is an important and not yet concluded research line. In this sense, it is important to note a significant increase in recent years in the number of remote sensing image segmentation-related scientific publications, together with the available image segmentation methods (e.g., semantic segmentation) [32]. Perhaps the introduction of SWIR bands in the MRS should also be studied in the future.
Regarding the classification accuracy assessments, both object-based (Table 3) and pixel-based (Table 4) tests reached similar results for all the three strategies, though the final pixel-based classifications were always slightly better. In previous works, the accuracies dealing with a sample of meaningful objects were always better than the pixel-based results achieved in the whole study area using a raster ground truth [21,22]. In the current work, and unlike the works cited above, the sample of objects used in the object-based accuracy assessments (i.e., 2350) covered most of the existing objects in the segmentation of the study area (about 3188). Moreover, as the segmentation of PCG was based on the ground truth, many of the misclassification problems due to errors in the segmentation were avoided.
The OA values significantly improved from 90.85% (using VNIR strategy) to 97.38% (with All Features strategy). The last OA value achieved by using the information contained in the 16 WV3 bands (VNIR plus SWIR) turns out to be the best that has been published in recent years for PCG or PMF [22,33,34,35,36]. It is important to note here that an even better OA value of 98.08% was reached using a knowledge-based DT for dealing with the whitewashed issue related to the disturbance of the spectral reflectance characteristics of the plastic films (Figure 7). With this strategy, only two PCG were wrongly classified as Non-GH. One of these greenhouses was made of glass, as already mentioned, while the other was a PCG abnormally whitewashed with a massive amount of lime.
It is important to note that Levin et al. [17] stated that the three main absorption wavelengths for agricultural plastic films were located at 1218 nm, 1732 nm and 2313 nm, and they were not affected by settling dust, whitewashing or by the underlying surface. Looking at Figure 7, the valley values reported by [17] could be appreciated only in the greenhouses without whitewashing.
The WV3 SWIR bands information, particularly in the form of NDPI and PMLI indices, provided valuable information to improve the separability between GH and Non-GH classes. However, the mean and SD object values of SWIR bands had a low importance in the PCG classification (Table 4). Furthermore, RBD index was not between the most important features for mapping PCG. In fact, any relationship between SWIR11 and SWIR13 bands can be found in Figure 7.
The result achieved in the study area was good thanks to the eight WV3 SWIR bands and the plastic discrimination power of NDPI. Unfortunately, WV3 is a commercial VHR satellite which is not available for free; this is a particular problem when applied to very large areas. In that sense, it would be necessary to study the capability of other optical sensors containing SWIR bands, such as Sentinel-2 or Landsat 8, or even the hyperspectral Precursore IperSpettrale della Missione Applicativa (PRISMA) sensor developed by the Italian Space Agency [37] for mapping PCG while taking into account the special characteristics of plastic materials.
Having verified the importance of the WV3 SWIR bands in the detection of PCG, we propose to go a step further, trying to distinguish between different types of plastic films used in the different PCG. According to [1], the main plastic films used in Spain are polyethylene (PE), polyvinylchloride (PVC) and polypropylene (PP), in that order. Moreover, up to 10% of greenhouse films are made of ethylene-vinyl acetate (EVA), sometimes being co-extruded with three layers made of EVA and low-density polyethylene (LDPE).
7. Conclusions
This study aimed to test the capability of the eight SWIR bands included in WV3 imagery (3.2 m GSD) to improve classification accuracy for mapping PCG attained from the VNIR eight bands (1.2 m GSD). For that, three strategies (VNIR, SWIR and All Features) were carried out by applying an OBIA approach.
The accuracy results clearly showed that the information available from WV3 SWIR bands (mainly NDPI and, secondly, PMLI indices) significantly improved the OA and kappa measures. Values of 90.85% and 0.812 for OA and kappa, respectively, were achieved using VNIR strategy, while improvements of around 6.7% for OA and 14% for kappa were reached working with the All Features strategy.
The NDPI index showed its potential for detecting plastic materials, resulting in the most important feature for mapping PCG. However, the fact that some PCGs within the study area were whitewashed in June 2020 decreased the NDPI values by masking the presence of plastic films. The last issue was addressed in this work by proposing a new index named NDPI_B that combines the NDPI index and Brightness for improving the detection of whitewashed PCG. Impressive accuracy values of 98.08% and 0.959 for OA and kappa, respectively, were obtained from using only the NDPI_B index.
Although the results achieved here are promising, we have only scratched the surface of the possibilities that WV3 offers for mapping PCG. Further works should be carried out in this research area.
Author Contributions
M.A.A. proposed the methodology, designed and conducted the experiments and contributed extensively to manuscript writing and revision. R.J.-L. collected the input data and performed the experiments, additionally collaborating on manuscript revision. F.J.A. made significant contributions to the research design and manuscript preparation. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by Spanish Ministry for Science, Innovation and Universities (Spain) and the European Union (European Regional Development Fund, ERDF) funds, grant number RTI2018-095403-B-I00.
Acknowledgments
We acknowledge the support of the Spanish Ministry for Science, Innovation and Universities (Spain) and the European Union (European Regional Development Fund, ERDF) funds (Grant Reference RTI2018-095403-B-I00). This study takes part in the general research lines promoted by the Agrifood Campus of International Excellence ceiA3 (
Conflicts of Interest
The authors declare no conflict of interest.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figures and Tables
Enlarge this image.Figure 1. Location of the study area: (a) General situation in Spain; (b) RGB WorldView-3 Multispectral orthoimage; (c) ground truth including greenhouse class (GH) in red and non-greenhouse class (Non-GH) in green. Coordinate system: WGS84 UTM Zone 30N.
Enlarge this image.Figure 3. Segmentation stage on a detailed area: (a) Ideal MRS segmentation achieved using the modified Euclidean Distance 2 (ED2); (b) manual digitized ground truth; (c) final segmentation including GH ground truth geometry.
Enlarge this image.Figure 4. Decision tree (DT) model proposed for All Features strategy (greenhouse class (GH), non-greenhouse class (Non-GH)).
Enlarge this image.Figure 5. NDPI index performance: (a) RGB WV3 orthoimage with GH polygons in red (the letters B in blue indicate agricultural buildings); (b) NDPI index and ellipses point out plastic areas belonging to Non-GH class; (c) green ellipse showing a plastic farm building; (d) blue ellipse depicting a plastics recycling plant; (e) white and yellow ellipses illustrate piles of plastic crates and plastic solar panels, respectively; (f) orange ellipse showing an abandoned greenhouse. All the detailed images come from Google Earth.
Enlarge this image.Figure 6. Pixel-based accuracy assessment over a detailed area (2000 m × 2000 m): (a) Classification results using the DT model for All Features strategy; (b) classification results using the NDPI_B index. The GH class is represented in red and the Non-GH class is represented in green. Yellow (GH classified as Non-GH) and purple (Non-GH classified as GH) colors in both classification maps highlight the misclassification errors.
Enlarge this image.Figure 7. Reflectance values of the SWIR bands of the WV3 image for four PCGs with different colors. The four orange vertical lines represent the central wavelengths of SWIR10, SWIR12, SWIR13 and SWIR16.
WorldView-3 PAN, VNIR and SWIR bands and their wavelength ranges [5].
| Band | Band Name | Image | Wavelengths (nm) | Band | Band Name | Image | Wavelengths (nm) |
|---|---|---|---|---|---|---|---|
| PAN | PAN | 450–800 | |||||
| 1 | Coastal (C) | VNIR | 400–450 | 9 | SWIR9 | SWIR | 1195–1225 |
| 2 | Blue (B) | 450–510 | 10 | SWIR10 | 1550–1590 | ||
| 3 | Green (G) | 510–580 | 11 | SWIR11 | 1640–1680 | ||
| 4 | Yellow (Y) | 585–625 | 12 | SWIR12 | 1710–1750 | ||
| 5 | Red (R) | 630–690 | 13 | SWIR13 | 2145–2185 | ||
| 6 | Red Edge (RE) | 705–745 | 14 | SWIR14 | 2185–2225 | ||
| 7 | NIR1 | 770–895 | 15 | SWIR15 | 2235–2285 | ||
| 8 | NIR2 | 860–1040 | 16 | SWIR16 | 2295–2365 |
WorldView-3 VNIR and SWIR object-based features. The WV3 band names (as given in Table 1) represent the related ground reflectance values.
| Source | Name | Definitions | References |
|---|---|---|---|
| VNIR | Mean and standard deviation (SD) (16 features) | Mean and SD of each of 1 to 8 band | [26] |
| Brightness | Average value for all 8 VNIR bands | [26] | |
| NDVI (Normalized Difference Vegetation Index) | (NIR2 − R)/(NIR2 + R) | [27] | |
| PGI (Plastic Greenhouse Index) | 100 × (B × (NIR2 − R))/(1 − (B + G + NIR2)/3) | [11] | |
| PI (Plastic Index) | NIR2/(NIR2 + R) | [9] | |
| SWIR | Mean and standard deviation (SD) (16 features) | Mean and SD of each of 9 to 16 band | [26] |
| NDPI (Normalized Difference Plastic Index) | ((SWIR10 − SWIR12) + (SWIR13 − SWIR16))/(SWIR10 + SWIR12 + SWIR13 + SWIR16) | [14] | |
| PMLI (Plastic–Mulched Landcover Index) | (SWIR10 – R)/(SWIR10 + R) | [12] | |
| RBD (Relative-absorption Band Depth index) | (SWIR11 − SWIR13)/SWIR13) | [15] |
Object-based confusion matrices for the three tested strategies.
| Strategy | User/reference | GH | Non-GH | Sum |
|---|---|---|---|---|
| All Features | GH | 1147 | 28 | 1175 |
| Non-GH | 51 | 1124 | 1175 | |
| Sum | 1198 | 1152 | 2350 | |
| OA (%) = 96.64 | kappa = 0.933 | |||
| VNIR | GH | 1034 | 141 | 1175 |
| Non-GH | 103 | 1072 | 1175 | |
| Sum | 1137 | 1213 | 2350 | |
| OA (%) = 89.62 | kappa = 0.792 | |||
| SWIR | GH | 1138 | 37 | 1175 |
| Non-GH | 79 | 1096 | 1175 | |
| Sum | 1217 | 1133 | 2350 | |
| OA (%) = 95.06 | kappa = 0.901 | |||
The top 10 object features for the three strategies, based on their relative importance.
| All Features | VNIR | SWIR | |||
|---|---|---|---|---|---|
| Feature | Importance | Feature | Importance | Feature | Importance |
| NDPI | 1.000 | Mean C | 1.000 | NDPI | 1.000 |
| PMLI | 0.968 | Mean B | 0.926 | PMLI | 0.956 |
| Mean C | 0.818 | PI | 0.862 | SD SWIR12 | 0.439 |
| PGI | 0.779 | NDVI | 0.862 | SD SWIR11 | 0.425 |
| Mean B | 0.759 | PGI | 0.781 | SD SWIR14 | 0.418 |
| Brightness | 0.711 | Mean G | 0.780 | SD SWIR10 | 0.417 |
| Mean G | 0.705 | Mean Y | 0.671 | SD SWIR13 | 0.413 |
| Mean RE | 0.660 | Brightness | 0.659 | SD SWIR15 | 0.411 |
| Mean Y | 0.647 | Mean R | 0.646 | SD SWIR16 | 0.403 |
| Mean R | 0.633 | SD R | 0.609 | Mean SWIR16 | 0.395 |
Pixel-based confusion matrices from ground truth (TTA Mask) for the three strategies.
| Strategy | User/reference | GH | Non-GH | Sum |
|---|---|---|---|---|
| All Features | GH | 8,691,478 | 97,357 | 8,788,835 |
| Non-GH | 278,134 | 5,271,465 | 5,549,599 | |
| Sum | 8,969,612 | 5,368,822 | 14,338,434 | |
| OA (%) = 97.38 | kappa = 0.944 | |||
| GH UA (%) = 98.89 | GH PA (%) = 96.90 | |||
| VNIR | GH | 7,791,948 | 134,237 | 7,926,185 |
| Non-GH | 1,177,664 | 5,234,585 | 6,412,249 | |
| Sum | 8,969,612 | 5,368,822 | 14,338,434 | |
| OA (%) = 90.85 | kappa = 0.812 | |||
| GH UA (%) = 98.31 | GH PA (%) = 86.87 | |||
| SWIR | GH | 8,622,938 | 113,534 | 8,736,472 |
| Non-GH | 346,674 | 5,255,288 | 5,601,962 | |
| Sum | 8,969,612 | 5,368,822 | 14,338,434 | |
| OA (%) = 96.79 | kappa = 0.932 | |||
| GH UA (%) = 98.70 | GH PA (%) = 96.14 | |||
Pixel-based confusion matrices from ground truth (TTA Mask) for the DT based on the NDPI_B index.
| Strategy | User/reference | GH | Non-GH | Sum |
|---|---|---|---|---|
| NDPI_B | GH | 8,905,137 | 211,147 | 9,116,284 |
| Non-GH | 64,475 | 5,157,675 | 5,222,150 | |
| Sum | 8,969,612 | 5,368,822 | 14,338,434 | |
| OA (%) = 98.08 | kappa = 0.959 | |||
| GH UA (%) = 97.68 | GH PA (%) = 99.28 | |||
© 2021 by the authors.
Details
