1.

Introduction

Cover crops are an essential agricultural best management practice benefiting growers. Cover crops are becoming more widely used by farmers, improving the sustainability of agriculture in the US.¹ Cover crops can contribute to the reduction of soil erosion, phosphorus loss, and groundwater nitrate leachate risk during winter.^2–5 Also, cover crops can protect perennial crop seedlings during establishment,^6,7 add organic matter, and enhance soil aggregation.⁷ In addition, leguminous cover crops can transform atmospheric nitrogen into biomass, which mineralizes in soil for the following grain crop to use.⁷ However, understanding the adoption and impacts of cover crops is hindered by the difficulty of accurately characterizing the spatial distribution of cover crops.

According to Ref. 8, cover crops are generally defined as grasses, legumes, and forbs planted for seasonal cover and associated benefits. These crops are meant to cover and enrich the soil instead of being harvested. The most common types of cover crops are rye and winter wheat.⁸ For example, cereal rye is often planted as a winter cover crop in the fall between cash crops (e.g., corn and soybeans).⁸ However, we focus on wintertime green vegetated cover (WGVC), which is overall vegetation cover on croplands during wintertime. The reason is that they function similarly to cover crop regardless of the purpose or precise timing of planting even though some of WGVC is not conventionally defined as cover crop. Any vegetation during wintertime is expected to have the ecological functionality of cover crops to a certain degree (such as reducing soil erosion), and perennial crops have been included in types of cover crop in some past studies (e.g., Refs. 9–11). Consequently, WGVC identified in this study may include perennial crops (e.g., alfalfa), field crops planted in fall or early spring (before May), and green weeds in addition to conventional cover crops. Therefore, we expected to have a larger area of WGVC compared to the surveyed area of cover crops of traditional definition.

Currently, the usage of cover crops is examined by windshield surveys^5,12 or by questionnaire surveys, which are incomplete, often biased, and require significant time, labor, and cost.¹³ Therefore, remote sensing technologies are known to be promising approaches for the assessment of cover crops especially in large scale cropland as noted in Ref. 13. A pixel-based approach has been applied in many cases (e.g., Refs. 5 and 13–15) to estimate cover crop area in agricultural remote sensing studies. Pixel-based approaches were used with agricultural and vegetation indices, such as a combination of normalized difference vegetation index (NDVI), normalized difference tillage index (NDTI), and normalized difference residue index (NDRI),⁵ or a combination of NDVI, difference vegetation index (DVI), normalized green red difference index (NGRDI), and ratio vegetation index (RVI).¹³ However, it is known that the pixel-based approach has a salt-and-pepper effect,^16–19 which can cause problems when converting classified pixels to polygons. In Ref. 13, it was noted that the salt-and-pepper effect was observed in pixel-based classification for the assessment of cover crops. It is obvious that this salt-and-pepper effect will occur in WGVC as well. Therefore, another remote sensing approach other than the pixel-based method was sought to overcome the salt-and-pepper effect when assessing WGVC. Object-based classification was applied in this study to detect WGVC from satellite imagery while solving the salt and pepper effect.

In many cases, object-based classification is applied using commercial fee-based software. However, an object-based classification workflow using commercial software financially burdens many environmental, agricultural, or educational organizations because of the high cost of software licenses. Therefore, the large-scale area WGVC detection workflow proposed in this study was implemented using open-source remote sensing and geographic information system (GIS) technologies. This study is feasible because of the growing maturity of open-source geospatial technologies. Related to this study, there has been an unmanned aerial vehicle-based object-based image analysis (OBIA) of cover crop detection in vineyards²⁰ using eCognition Developer 9.2²¹ recently. However, a workflow for large-scale (e.g., state-level) WGVC identification by supervised object-based classification using open-source technologies was not found in the literature.

For large-scale (e.g., state-level) cropland WGVC detection, it is ideal to use medium resolution imagery. Processing a high-resolution image for a large area will require large computing resources, and low-resolution images will give inaccurate results for WGVC delineation. Sentinel-2 constellation images have 10 m resolution with 5-day temporal resolution. Also, these images are free of charge to public, scientific, and commercial users.²² Therefore, Sentinel-2 images were used in this study because they are mid-resolution satellite images with near real-time accessibility.

The purpose of this study is to develop a workflow for delineating large-scale WGVC with OBIA and open-source geospatial technologies using Delaware, United States as our area of study. Delaware is a good location because it allows us to process and analyze data for an entire small state with significant agricultural area. Hence, this study will fill the absence of OBIA studies for large scale identification of WGVC using open-source geospatial technologies.

2.

Materials and Methods

The selected area of interest is Delaware, United States. For this area of interest, Copernicus Sentinel-2^23,24 (specifications in Table 1) images were acquired from US Geological Survey (USGS) EarthExplorer.²⁵ Currently, the Sentinel-2 data are accessible in Copernicus Browser²⁴ but not in USGS EarthExplorer any longer. The revisit time of 5 days was sufficient to find cloudless images in reasonable temporal intervals.

Table 1

Description of Sentinel-2 constellation.23,24

Because WGVC can change during a season, multiple image dates were needed to represent the entire wintertime. The selected images were from December of the harvest year, and from February and April of the following year. Since three tiles are required to cover Delaware, virtual raster images were built using the raster program gdalbuildvrt²⁶ using tiles in Table 2 and clipped in “.tif” format for the state of Delaware.

Table 2

Description of image data.

These false color (NIR, R, and G bands) images (Fig. 1 displayed with cumulative count cut with minimum of 2% and maximum of 98%) were segmented using large-scale mean-shift (LSMS)²⁷ segmentation with Orfeo ToolBox (OTB).^28,29 The mean shift algorithm³⁰ has advantages in its use of a hierarchical relationship between segmentation levels, unlike the scale invariance of the watershed algorithm.³¹ Also, it does not require prior knowledge about cluster numbers and shape constraints.³² In addition, more complex shapes can be modeled using the mean shift algorithm compared to K-means.³³ LSMS was chosen as the segmentation algorithm because the image data in this study had 10 m GSD and covered the entire state of Delaware, which was quite large to be segmented using the traditional mean shift option in OTB. Table 3 shows the parameters used for LSMS in this study.

Fig. 1

Sentinel-2 false color images (NIR, R, and G) of 12/26/2021, 2/9/2022, and 4/30/2022 (from left to right), shown using cumulative count cut with a minimum (2%) and maximum (98%).

Table 3

Parameters for large scale mean shift segmentation.

Each spatial and range radius in Table 3 represents the thresholds of spatial and spectral signature Euclidean distance in evaluating pixels in an object.³⁴ The minimum size for segmentation of an object was set to 400 pixels, which is equivalent to 9.88 acres. This value was determined after finding that the size of most individual fields was larger than 400 pixels after visual examination of the fields in Sentinel-2 images. Also, any isolated crop field less than 400 pixels would be defined by later refinement—extraction of cropland area using 2021 National Agricultural Statistics Service (NASS) cropland data layer (CDL)³⁵ boundary.

To restrict the area of interest to field crop areas, only the segmentations on the field crops were extracted by clipping segmented objects by CDL cropland areas. Beyond a small-scale area (e.g., farm level), it is hard to identify field crops physically for the entire state. Therefore, 2021 NASS CDL for Delaware was used to identify field crops area, as used for cover crop identification by Hively et al.^5,36 in a similar way. Because CDL includes land covers other than crops, only field crop pixels were extracted using gdal_calc.py. Table 4 shows the field crops and their area calculated using QGIS³⁷ raster layer unique values report tool. The crop name for each pixel value could be found in Delaware CDL Metadata.³⁸ The total area of field crops was calculated as 183,754 ha (454,057 acres). The most significant field crops were corn and soybeans, whose sum represented about 81% of the entire field crop area.

Table 4

Candidate summer crop area from 2021 Delaware CDL for WGVC. (Abbreviations: Dbl = double, WinWht = winter wheat.)

However, the clipping by CDL modified the segmentation object polygons and even created multi-part polygons in some locations after clipping out central areas of polygons. Therefore, the attribute of each object had to be recalculated to reflect these changes. First, possible multipart polygons were separated by using Multipart to Single part tool in QGIS, and the objects’ attributes (mean and variance) were calculated with the Zonal Statistics tool in QGIS for each band. In Ref. 39, texture was measured in terms of homogeneity, contrast, dissimilarity, entropy, standard deviation, correlation, angular second moment, and mean. Therefore, using mean and variance addresses not only spectral characteristics but also textures of objects to some degree.

After training samples, the field crop segmentations were classified using normal Bayes (NB)⁴⁰ supervised classification in OTB. In Ref. 41, NB and support vector machine (SVM) showed better performance than the classification and regression tree (CART) and K nearest neighbor (KNN). Although NB needs a higher number of samples, SVM requires complex tuning parameters. Therefore, NB has a clear advantage in this study compared to SVM, CART, and KNN. For sampling purposes, both false and true color images, which were displayed with cumulative count cut with minimum (2%) and maximum (98%), were used as ground truth. Using best judgment, the operator could identify most vegetation cover on the field crop area as red and green color tones in each false and true color image for training. Figure 2 shows examples of the sampling of WGVC and non-WGVC on the clipped segmentations for training. The sampling points for each class were marked inside each sample object by the operator. The properties of each object were acquired for each sample point by join attributes by location with intersects geometric predicate in QGIS.

Fig. 2

Sampling example of WGVC (green dot) and non-WGVC (red dot).

It was noted that the number of training samples should be between 10 and 100 times the number of bands in practice.⁴² Therefore, the number of samples, which gave satisfactory classification results later in each case, were larger than 30 (Table 5) in every case. This operation was implemented for each false color image for each of the three dates, and they were combined by Union operation in QGIS to represent the WGVC between the winter 2021 and the spring of 2022.

Table 5

The number of training samples required for classification results to become satisfactory for each image date.

The average NDVI values of WGVC and non-WGVC objects for each date were compared. After creating NDVI raster, the mean NDVI for each sample object was calculated by zonal statistics in QGIS.³⁷ In Fig. 3, WGVC samples had higher NDVI values than non-WGVC samples. However, the minimum NDVI of WGVC and maximum NDVI of non-WGVC overlapped to some degree for 12/26/21 and 02/09/22. The NB classifier is expected to address this confusion problem by computing probability of membership to each class compared to alternative methods relying on NDVI thresholds.

Fig. 3

Boxplot of NDVI of WGVC and non-WGVC samples from (a) 12/26/21, (b) 02/09/22, and (c) 04/30/22. Boxes show interquartile range with solid lines for medians and dotted lines for means, and whiskers show minima and maxima (or upper fences when a point represents a maximum).

Figure 4 shows the main workflow described above.

Fig. 4

WGVC area estimation workflow.

For the accuracy assessment, a different dataset other than Sentinel-2 imagery was sought. The criteria of the data for accuracy assessment were (1) the data should be publicly available geospatial data; (2) the data should show land cover on the dates near those of Sentinel-2 images. Landsat 8 and 9 images satisfied the criteria, since they are publicly available georeferenced images in USGS EarthExplorer, and Landsat has temporal resolution frequent enough to have similar image capture dates with Sentinel-2 images. Using Landsat 8 and 9 images, false color images were composited with NIR, R, and G bands, and true color images were composited with R, G, and B bands, as shown in Table 6. Both false and true color images were used to assess the accuracy of the classification. However, the GSD of NIR, R, G, and B bands of Landsat 8 and 9 were quite large (30 m). Therefore, pan-sharpened images ($\mathrm{GSD}=15\mathrm{m}$) were created for both false and true color composite images using panchromatic images ($\mathrm{GSD}=15\mathrm{m}$) with Pansharpening tool⁴³ in the GDAL plugin of QGIS. Figure 5 shows the images used for the accuracy assessment. The time difference between Landsat and Sentinel data was 10 days maximum. This time gap was allowed because some of the Landsat 8 and 9 imagery taken near Sentinel-2 data acquisition dates had severe amounts of clouds. Fifty samples with an area larger than or equal to 10,000 m² were collected randomly from each class of classification results using QGIS random extract within subsets tool for date 12/26/2021 and 2/9/2022 results. The classification results of those objects were compared with pan-sharpened false color and true color Landsat images (12/26/2021, 2/9/2022 each) by the operator.

Table 6

Landsat data used for accuracy assessment.

Fig. 5

(a) False color Landsat 8 (12/16/21), (b) false color Landsat 9 (2/10/22), (c) false color Landsat 8 (5/9/22), (d) true color Landsat 8 (12/16/21), (e) true color Landsat 9 (2/10/22), and (f) true color Landsat 8 (5/9/22).

However, the 5/9/2022 Landsat 8 image still had thick cloud cover in the southeastern part of DE. Therefore, 80 initial sample polygons larger than 10,000 m² were randomly chosen from 4/30/2022 classification results table, and 50 samples, which were not affected by the cloud cover in the order of rows were used for the accuracy assessment by comparing with pan-sharpened false and true color Landsat image of 5/9/2022. The area threshold (10,000 m²) was imposed to include only larger polygons for accuracy assessment. If the sample object was composed of WGVC and non-WGVC pixels, the object was identified as the class of the majority of pixels. Since the samples for the accuracy assessment were segmented objects, the classification accuracy was calculated considering the area of each sample object using Eq. (1) according to Ref. 44 (as cited in Ref. 45):

Also, the Kappa statistic was calculated with the area of objects instead of the number of pixels in the following equation:⁴²

The strength of the Kappa statistics was evaluated for agreement as suggested in Ref. 46, with values of 0.41 to 0.60 considered moderate, 0.61 to 0.80 considered substantial, and 0.81 to 1.00 considered almost perfect agreement. Finally, the areas of WGVC for each date and the area merged into polygons by Union, which is a QGIS vector processing tool, were calculated.

3.

Results and Discussion

OTB provided the results of training of samples in the form of a confusion matrix (Table 7) for WGVC and non-WGVC. The training accuracy of every image was high.

Table 7

Confusion matrix for training sampling (column is reference label, row is produced label).

WGVC and non-WGVC of Delaware were classified as shown in Figs. 6(a)–6(c) for 12/26/2021, 02/09/2022, and 04/30/2022 using an NB classifier with trained models for the extracted objects for field crop area. Also, visual inspection suggested WGVC and non-WGVC usually was classified properly (example zoomed area west of Dover, Delaware, United States in Fig. 7). In Fig. 6, there was no substantial WGVC in the northern part of DE, because this part is mostly urban area (e.g., the cities of Newark and Wilmington). Also, the Delaware Bay area (mostly wetlands), Dover (urban), and Redden State Forest (forest) did not have a substantial presence of crops. The polygons merged by union for the above three dates showed WGVC identification spread across most of the state [Fig. 6(d)]. The accuracy of classification was evaluated against pan-sharpened Landsat 8, 9 false/true color images by visual inspection as shown in the confusion matrices (Table 8).

Fig. 6

WGVC classification: (a) classification - 12/26/2021, (b) classification - 02/09/2022, (c) classification - 04/30//2022, and (d) WGCV merged by Union.

Fig. 7

False color image (a) and classification (b) for an example location west of Dover, Delaware, United States.

Table 8

Confusion matrices and Kappa statistics of three classification result.

In confusion matrices, user’s, producer’s, and overall accuracies in all dates were higher than 85% (Table 8). Therefore, the classification results were entirely satisfactory. Also, kappa statistics (KHAT), which is a measure of the difference between the actual agreement and chance agreement,⁴² was higher than 74% in all cases. These kappa statistics were quite satisfactory, as well. Total WGVC area was found to be about 137,297 ha (Table 9).

Table 9

The total area of WGVC class union.

It was found that a multi-date WGVC identification strategy was useful because the total area estimated by Union was higher than found using images of individual dates. The total area of field crops from Delaware CDL 2021 (Table 4) was 183,754 ha. Therefore, about three quarters (75%) of the field crop areas were covered by WGVC between 12/26/2021 and 4/30/2022. Since WGVC includes conventionally defined cover crop, it is meaningful to compare WGVC area of 2021 to 2022 winter with the cover crop area most recently surveyed in 2022. In the Nonpoint Source Program (NPSP) 2022 Annual Report for Delaware,¹² the cover crop area was 40,811 ha (originally 100,846 acres) or $\sim 22\%$ of field crops. When we compare this with total WGVC found (137,297 ha), the difference is 96,486 ha. The first reason of difference is the additional types of vegetation of WGVC that are not defined as conventional cover crop. The second reason is the difference in traditionally defined cover crop area itself between the survey and true value, since the survey methods can be incomplete and biased.¹³

It should be noted that the usage of cover crops has been found to be increased during the investigation. The cover crop usage has been increased by 50% between 2012 and 2017 in the United States,⁸ and four times between 2011 and 2021 in the Midwestern United States.⁴⁷ The 2016 to 2017 report on the fifth annual cover crop survey by the Sustainable Agriculture Research and Education and the Conservation Technology Information Center, where respondents represented 47 states, showed about 25% and 60% increase in the usage of cover crops from 2014 to 2015 and 2016, respectively.⁴⁸ Also, the 2022-2023 report, where respondents represented 49 states, showed that the mean number of acres of cover crop for respondents who used cover crops, has been increased from 324.9 acres to 413.6 acres between 2018 and 2022.⁴⁹ This upward trend was verified when comparing the cover crop area surveyed in the NPSP 2022 Annual Report for Delaware¹² and the Census of Agriculture (COA) of 2017⁵⁰ for Delaware. The cover crop area of Delaware in the 2022 NPSP annual report was 40,811 ha, which was larger by 5,153 ha when compared with the cropland planted to a cover crop (excluding CRP) (35,658 ha, originally 88,112 acres) in 2017 COA for Delaware. Because vegetation during winter included not only planted cover crops, the WGVC estimated by the proposed method in this study would be more beneficial for some agricultural and environmental modeling applications.

Using OTB for supervised OBIA is challenging because of the lack of user-friendly documents for sampling, training, classification, and accuracy assessment. This is often the case for many free, open-source technologies. However, the object-based classification tutorial⁵¹ was informative for forming the workflow of supervised classification using OTB in this study. The current study is expected to contribute to these resources, since we conducted supervised object-based classification with concrete examples in this study. Therefore, this study is applicable to workflow development for other uses of OBIA.

4.

Summary and Conclusion

To identify WGVC areas in Delaware, object-based classification was applied using open-source geospatial technologies. The application of a remote sensing technique (object-based image classification) with LSMS segmentation enabled large-scale WGVC detection, which is efficient in terms of cost, time, and labor. Also, the salt-and-pepper effect was removed by applying an object-based classification approach instead of the traditional pixel-based methodology. Another hurdle in object-based classification, which is fee-based commercial software’s high cost, was overcome using OTB, GDAL, and QGIS, which are free and open-source geospatial technologies. For training NB classifier, samples were chosen for WGVC and non-WGVC classes. For accuracy assessment, classification results were compared with the pan-sharpened Landsat 8,9 images by the operator. The final NB classification results evaluated against pan-sharpened Landsat 8,9 false/true color images were quite satisfactory with confusion matrices showing overall accuracies higher than 85% and KHATs higher than 74% in all cases. The input data was created using NIR, R, and G bands of Sentinel-2 images. These images were free and publicly available with optimal GSD (10 m) for state-level analysis and high temporal resolution (5 days), which enabled the acquisition of cloudless images in a reasonable time interval. The field crop areas from NASS CDL were used as candidate areas of WGVC. The total WGVC area was created by using the Union tool in QGIS, merging cover crop polygons of three dates (12/21/2021, 02/06/2022, and 04/30/2022), and this area was estimated as 137,297 ha overall.

The presented study offers large-scale (state-level) WGVC detection with supervised object-based classification using completely open-source technologies. The presented workflow in this study will be beneficial to future WGVC studies and can be used by organizations with limited time, labor, and funding.