Crop-Type Classification for Long-Term Modeling: An Integrated Remote Sensing and Machine Learning Approach

Long-term temporal and spatial information of crop type supports a wide range of applications including hydrological and climatological studies. In the U.S., yearly crop data layers (CDLs) are available starting in the early 2000s and have been developed using combined field information and sets of temporal imagery from multiple sensors. Development of long-term crop-type layers similar to CDLs is restricted by reduced accessibility to imagery and the necessary auxiliary datasets. In this study, a procedure to generate a historical crop type was developed and evaluated. Time series of Normalized Difference Vegetation Index (NDVI) datasets from Landsat 5 TM sensor for the Lower Bear Creek watershed were collected and processed. Object-based pseudo phenology curves, represented by the NDVI time series, were generated using noise filtering and dimensionality standardization procedures for the years 1985, 1990, 1995, 2000, and 2005. Classifiers were developed and evaluated using random-forest machine learning algorithms and CDL datasets as the reference. Increased generalization performance was obtained when the model was developed using multi-year datasets. This can be attributed to improved crop type representation during the training phase coupled with characterization of yearly variations due to natural (weather) and anthropogenic factors (farming management). Source of uncertainties were the presence of multiple crops within objects, phenological similarities between soybean and corn/maize, and the accuracy of CDL itself. The proposed procedure supports the development of historic crop types for long-term studies at the field scale in agricultural watersheds.


Introduction
Datasets on crop type classification and its distribution in space and time support a wide range of direct applications and can be incorporated into a variety of environmental and hydrologic models. The latter allow us to assess the impact of agriculture on the environment, and conversely, the impact of the environment on agriculture [1,2]. These datasets support government agencies, private sector organizations, scientists, educators, and others seeking an improved understanding of cropland spatiotemporal variations and make informed decisions on farming and conservation practices [2,3]. Applications include, but are not limited to, crop yield estimation and prediction [1], natural hazard (e.g., floods and droughts) impact assessment on agricultural commodities [4][5][6], and nonpoint source pollution quantification and mitigation [5,6]. For the latter, hydrological watershed models, in particular the annualized agricultural non-point source (AnnAGNPS) [7,8] and the Soil and Water Assessment Tool (SWAT) [9,10], rely heavily on spatiotemporal crop type information to characterize agricultural watershed and best represent existing conditions of land use and farming practices.
The commonality between these studies is their focus on developing detailed and accurate methodologies for crop type mapping for recent years (2000 to the present). Limited studies were found addressing historical crop type classification at the field scale based solely on remote sensing products. In this manuscript, we developed and evaluated a field scale crop type classification methodology that relies on limited remote sensing datasets (i.e., Landsat products) without requiring any auxiliary data. Limited temporal resolution is addressed by an integrated data processing and filtering method with robust machine learning algorithms. Hence, this methodology could readily be applied in many parts of the world as well as in the U.S. to complement (i.e., 1985-2013) the more recent timeframe covered by CDLs (i.e., [2006][2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017][2018][2019]. This could provide a historical distribution of cropland for a better understanding of long-term trends in agricultural commodities, yield and variations and their relation with climate change and other long term anthropogenic and natural controlling variables.

Study Site
The Lower Bear Creek watershed with an area of 93,833 hectares is located in the southeastern corner of the State of North Dakota, USA ( Figure 1). It is characterized by low relief (396 to 894 m a.m.s.l.) and land use/land cover predominantly of agriculture in the form of row crops and grazing fields (> 90%). According to the USDA cropland data layers, in 1998, the main composition land use/land cover was 32% fallow, idle, grass, or pasture, and 60% row crops. The main row crops in 1998 were corn/maize (10%), soybeans (11%), and spring wheat (22%). In 2005, these percentages changed to 14%, 31%, and 16% for corn/maize, soybeans, and spring wheat, respectively. Crop conversion changes aligned well with the national trends of conversion of traditional crops into soybeans and corn/maize. In North Dakota, typical planting and harvesting dates for corn/maize are 2-28 May and 8 October-19 November, for soybean 14 May-3 Jun3 and 24 September-21 October, and for spring wheat 24 April-25 May and 8 August-13 September [33].

Imagery Acquisition and Processing
Normalized difference vegetation index (NDVI) and pixel quality metadata (PQA) images were downloaded from the U.S. Geological Survey (USGS)-Earth Resource Observation and Science (EROS)-EROS Science Processing Architecture (ESPA) website [34,35]. All images in the years 1985, 1990, 1995, 2000, and 2005 were downloaded (Table 1). Only images from the Landsat 5TM collection were used because this mission covered all five years considered. The product Landsat surface reflectance-derived spectral indices (NDVI) was chosen to reduce pre-processing steps such as atmospheric correction and spectral index computation [35]. The study site is covered by two Landsat paths (Figure 1d). Both NDVI and PQA images were re-projected ( Figure 2, box 1) and clipped to the watershed extents ( Figure 2, box 2).
Parcel boundaries were digitized based on high-resolution aerial photographs and labeled as agriculture or non-agriculture, based on the visual inspection of aerial photographs and supporting auxiliary data such as CDL datasets for 2000 and 2005, road centerlines, and NLCD ( Figure 2, box 3). Out of 2757 parcels, 1755 were identified as croplands (Figure 1c). Cropland fields were used to generate a non-agriculture mask imagery that was subsequently applied to all NDVI images ( Figure 2, box 4). Similarly, a set of masks were developed to filter out the potential influence of clouds, cloud shadows, and snow/ice to the NDVI values using the flags encoded in the PQA images [36,37] and generated by the CFMask algorithm [35] (Figure 2, box 5 and Figure 3a). use/land cover was 32% fallow, idle, grass, or pasture, and 60% row crops. The main row crops in 1998 were corn/maize (10%), soybeans (11%), and spring wheat (22%). In 2005, these percentages changed to 14%, 31%, and 16% for corn/maize, soybeans, and spring wheat, respectively. Crop conversion changes aligned well with the national trends of conversion of traditional crops into soybeans and corn/maize. In North Dakota, typical planting and harvesting dates for corn/maize are 2-28 May and 8 October-19 November, for soybean 14 May-3 Jun3 and 24 September-21 October, and for spring wheat 24 April-25 May and 8 August-13 September [33].   Individual images of masked NDVI for cropland and cloud/snow were stacked into five files (one for each year considered). The final NDVI time series images contained a varying number of bands, each representing a day-of-year (DOY). Custom GIS zonal statistics algorithms were developed to calculate the average NDVI value for each DOY (each band) and each field (each image object). A modified zonal statistic algorithm was developed and employed to consider only valid pixels based on previously developed masks ( Figure 3a) and to report an average, if a minimum number of five or more valid pixels was found. This information was recorded in a tabular format of one field per row and the sequence of NDVI average values in columns ( Figure 2, box 6).
Reference datasets for 2000 and 2005 were derived from the respective CDLs. The raw CDL raster grid layers contained noise ( Figure 3b) and were processed through another custom-made zonal statistic algorithm. Common zonal statistics tools found in off-the-shelf GIS software packages calculate the results using all pixels and report only the class number with the majority value found. The algorithm employed recorded the frequency information of all classes found in each parcel, rather than only the most abundant one (Figure 3c). The dominant crop-type for each field was added to the NDVI temporal sequence database generated previously.    Individual images of masked NDVI for cropland and cloud/snow were stacked into five files (one for each year considered). The final NDVI time series images contained a varying number of bands, each representing a day-of-year (DOY). Custom GIS zonal statistics algorithms were developed to calculate the average NDVI value for each DOY (each band) and each field (each image object). A modified zonal statistic algorithm was developed and employed to consider only valid pixels based on previously developed masks ( Figure 3a) and to report an average, if a minimum number of five or more valid pixels was found. This information was recorded in a tabular format of one field per row and the sequence of NDVI average values in columns ( Figure 2, box 6).
Reference datasets for 2000 and 2005 were derived from the respective CDLs. The raw CDL raster grid layers contained noise ( Figure 3b) and were processed through another custom-made zonal statistic algorithm. Common zonal statistics tools found in off-the-shelf GIS software packages calculate the results using all pixels and report only the class number with the majority value found. The algorithm employed recorded the frequency information of all classes found in each parcel, rather than only the most abundant one (Figure 3c). The dominant crop-type for each field was added to the NDVI temporal sequence database generated previously. Individual images of masked NDVI for cropland and cloud/snow were stacked into five files (one for each year considered). The final NDVI time series images contained a varying number of bands, each representing a day-of-year (DOY). Custom GIS zonal statistics algorithms were developed to calculate the average NDVI value for each DOY (each band) and each field (each image object). A modified zonal statistic algorithm was developed and employed to consider only valid pixels based on previously developed masks ( Figure 3a) and to report an average, if a minimum number of five or more valid pixels was found. This information was recorded in a tabular format of one field per row and the sequence of NDVI average values in columns ( Figure 2, box 6).
Reference datasets for 2000 and 2005 were derived from the respective CDLs. The raw CDL raster grid layers contained noise ( Figure 3b) and were processed through another custom-made zonal statistic algorithm. Common zonal statistics tools found in off-the-shelf GIS software packages calculate the results using all pixels and report only the class number with the majority value found. The algorithm employed recorded the frequency information of all classes found in each parcel, rather than only the most abundant one (Figure 3c). The dominant crop-type for each field was added to the NDVI temporal sequence database generated previously.

Noise Filtering and Dimensionality Standardization
Clouds, cloud shadows, snow, and ice prevent satellite sensors from directly measuring vegetation growth by obscuring and contaminating areas of interest. Removing cloudy pixels creates troughs (instances of zero growth) in the growth curve that misrepresent the natural vegetation development cycle. Furthermore, the pixel quality assurance band used to mask out images is not perfect. Cloudy or shadowed pixels may persist even after masking with the PQA image. In order to best capture a reliable phenology curve from the aggregated data, noise in the phenology curve must be removed to address the substantial variability inherent to our dataset.
The procedure adopted was the weighted windowed linear regression [38]. Weighted windowed linear regression attempts to address variability by smoothing phenology curves to remove noise. A weight is assigned to each day of the year based on whether the point is a peak point, upward sloping, downward sloping, or a valley point where a peak point receives the highest weight and a valley point receives the lowest weight. A moving window is then used to successively fit a linear equation to each point in the field's weighted phenology curve. The resulting coefficients are used to calculate a series of weighted NDVI values for each day of year, which are then averaged where peak values are weighted heavier than valleys. The resulting product is a smoothed phenology curve where the growth obscured previously by gaps is restored ( Figure 4a). This procedure was applied to all fields and all five years.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 16 Figure 3. Development of input and reference datasets for individual agricultural fields. Red patches represent clouds used to mask/remove contaminate pixels during the field spatial aggregation of NDVI pixel values (a). The raw crop data layer (CDL) dataset (b) was processed to determine the majority crop type for each field based on the pixel-field zonal statistics (c). This operation was performed for the 2000 (shown) and 2005 datasets.

Noise Filtering and Dimensionality Standardization
Clouds, cloud shadows, snow, and ice prevent satellite sensors from directly measuring vegetation growth by obscuring and contaminating areas of interest. Removing cloudy pixels creates troughs (instances of zero growth) in the growth curve that misrepresent the natural vegetation development cycle. Furthermore, the pixel quality assurance band used to mask out images is not perfect. Cloudy or shadowed pixels may persist even after masking with the PQA image. In order to best capture a reliable phenology curve from the aggregated data, noise in the phenology curve must be removed to address the substantial variability inherent to our dataset.
The procedure adopted was the weighted windowed linear regression [38]. Weighted windowed linear regression attempts to address variability by smoothing phenology curves to remove noise. A weight is assigned to each day of the year based on whether the point is a peak point, upward sloping, downward sloping, or a valley point where a peak point receives the highest weight and a valley point receives the lowest weight. A moving window is then used to successively fit a linear equation to each point in the field's weighted phenology curve. The resulting coefficients are used to calculate a series of weighted NDVI values for each day of year, which are then averaged where peak values are weighted heavier than valleys. The resulting product is a smoothed phenology curve where the growth obscured previously by gaps is restored (Figure 4a). This procedure was applied to all fields and all five years. Although the smoothed curves are a better representation of field conditions in time, the curve's dimensionality varied per year and even per fields within the same year. Fields were described by different DOYs and varying number of samples. This was due to differences in the DOY when images were acquired (see Table 1) and potential gaps were present in one field, but not in another. This was addressed through the function fitting approach to existing points using a piecewise cubic polynomial method [39] implemented in the Scientific Computing library of the Python programming language [40,41]. The advantage of this method is the generation of a smooth curve that goes through the observed points. Once the coefficients are determined, the fitted equation can then be used to estimate values at set time intervals. Values were estimated for every 10 days from DOY 100 to DOY 330 (Figure 4b), assuring all curves (all fields and years) contained the same number of points; a requirement for machine learning algorithms. Using the CDL information assigned to each field, averaged curves for prevailing crop types were developed for 2000 ( Figure 4c) and 2005 (Figure 4d).

Models Development and Evaluation
The development/utilization of robust classification algorithms for crop-type classification has relied on machine learning algorithms (non-parametric) [28]. These algorithms are designed to address complex problems without being explicitly programmed to do so, but rather to derive solutions from recursive and iterative analysis of candidate solutions based on a user-provided training dataset [42]. Non-parametric algorithms for time-series analysis have significant advantages. These algorithms have been demonstrated to be distribution independent, work well with missing and/or noise data, have a higher accuracy than standard remote sensing classifiers, and are considered robust (capable of generalization). Multiple machine learning algorithms were evaluated using the WEKA Machine Learning Software package [43] with the random forest classifier [44] selected due to its higher accuracy over the artificial neural network, decision trees, and support vector machine. Random forest is described as an ensemble data mining/classification algorithm because it evolves and evaluates a set of independent tree-based classifiers concurrently. During the development of candidate solutions, the algorithm relies on bootstrapping techniques for improved conversion by assuring diversity [32]. For classification problems, assignment of an instance to class is based on predictions from all trees (majority voting). The limitations of this approach is that the individual classification trees are not available to the user and it can be computationally intense. Conversely, this method has the advantages of being able to quantify the importance of each input variable, it yields higher accuracy than traditional decision trees, and can generate accurate robust classifiers [42].
In this study, each crop field (image object) for a specific year corresponded to an instance. Each instance was described by 23 processed NDVI values for pre-determined DOYs between DOY 100-330 (representing growing season) at a 10-day time step (input variable) plus label information, when available (for 2000 and 2005 only). A total of three models were developed using NDVI images and CDL from 2000, 2005, and combined 2000-2005. For models developed using a single year, training was performed using cross-validation with 10 folds and testing using the other year. For the combined 2000-2005 datasets, the dataset was randomly divided into 66% for training and 33% for testing. During the development, the random forest algorithm proposed by [44] and implemented in the WEKA software package [43] was used with the following controlling parameters: bag size of 100%, minimum numeric class variance proportion of train variance for split of 0.001, minimum number of variables per leaf of 1, and maximum number of iterations of 1000. Models were assessed using standard statistics of overall accuracy and kappa coefficient of agreement [45] derived from the confusion matrix.

Pseudo Crop Phenology Curves
Analysis of the average NDVI curves for each class and each time period (Figure 4c,d) indicates patterns within crops (same crop different time period) and between crops (different crops). There are phenological similarities between corn, soybeans, and sunflower crops. These temporal signature similarities could restrain the classifier's performance of discerning between them. The similarities were more pronounced in 2000 (Figure 4c Analysis of the average NDVI curves for each class and each time period (Figures 4c,d) indicates patterns within crops (same crop different time period) and between crops (different crops). There are phenological similarities between corn, soybeans, and sunflower crops. These temporal signature similarities could restrain the classifier's performance of discerning between them. The similarities were more pronounced in 2000 (Figure 4c   The lack of reference, prevented absolute accuracy assessment and therefore a relative accuracy assessment was performed by comparing the predictions between models. Evaluating the agreement amongst all three models for the predicted correctly planted area for all classes was 68.7%, 75.9%, and 65.46% for 1985, 1990, and 1995, respectively. Slightly higher agreement amongst the three models was observed for 1990. Considering the agreement between the pairs of models, higher agreement was obtained between models

Models Development and Evaluation
Single year models were comparable (Table 2) and yielded high overall accuracy and kappa coefficient of agreement values for both training and testing (kappa > 0.55 and OA > 0.67). However, despite the model's overall high performance, a misclassification of corn/maize, soybeans, and sunflower can be noted. This can be attributed to the phenological similarities between these crops (yellow circle, blue diamond, and green square curves in Figure 4c,d).

Model Application
The three models developed, 2000, 2005, and combined 2000-2005 were applied to process the NDVI time series for 1985, 1990, and 1995; years for which CDL datasets are not available ( Figure 6). The lack of reference, prevented absolute accuracy assessment and therefore a relative accuracy assessment was performed by comparing the predictions between models. Evaluating the agreement amongst all three models for the predicted correctly planted area for all classes was 68.7%, 75.9%, and 65.46% for 1985, 1990, and 1995, respectively. Slightly higher agreement amongst the three models was observed for 1990. Considering the agreement between the pairs of models, higher agreement was

Pseudo Crop Phenology Curves
Working in the object domain, rather than in the pixel domain, has advantages. First, working with image objects removes the need for geospatial registration (georectification) between images to assure that individual pixels align spatially (since temporal signatures are generated by stacking multiple scenes). Georectification is a process that is hard to automate, and therefore often relies heavily on manual labor that is time-consuming and can add subjectivity. Using an object-based approach for the temporal analysis of spectral indices provides spatially smoothed NDVI values (field mean), reducing the noise/influence from individual pixel outliers. It increases the availability of data when information from an object is partially compromised (field partially covered by clouds). Additionally, there is significant reduction in the data volume from going from pixel to object scales (from millions of pixels to hundreds of fields). Similar studies have found that working in the object domain rather than in the pixel domain reduces noise and, as result, yields improved accuracy results [27,46]. However, changing from pixel to object is not without its challenges. There is the need to define field boundaries. In this study, this was done by heads-up digitizing using CDL and aerial photographs from one time period, therefore, it is possible that the field boundaries changed over time. Furthermore, a field could contain more than one crop in a particular growing season, causing potential bias in the average NDVI values. This was addressed by changing from simple majority voting (standard GIS zonal statistics operation) to determining the percentage of each crop in each field, and then using only fields with a percent majority greater than 50% (custom-developed zonal statistic).
The final processed phenology curve used in the models did not contain gaps. These curves are generated by two procedures: noise filtering (weighted windowed linear regression) and dimensionality standardization (piecewise cubic polynomial). However, these datasets are not immune to noise. Inter-field variations could be potentially attributed to gaps in the dataset (clouds) incorrectly filled by subsequent steps, inconsistencies in the CDL layer used as the reference, or contaminated pixels (cloud shadows) affecting the individual time period. This uncertainty is often addressed with increased temporal resolution, but remains a challenge when the objective is to classify historic datasets where temporal resolution is limited.
Averaged curves generated for corn and soybeans were compared with state-wide five-year averaged curves generated using MODIS (Figure 7) [47]. For both crops in 2000, this comparison demonstrated a standard development stage, but with a shorter peak, followed by an earlier decline. The shaded area in Figure 7 describes the variance represented by one standard deviation. It can be observed a higher variance for 2000 between DOY 210 and DOY 260. The curves generated in 2005, despite a slightly late start, matched the state-wide averages. The late start can be attributed to abnormally dry weather between April 26, 2005 to May 31, 2005 [48] followed by high rainfall amounts in May/June. The comparisons of the developed with published curves generated using different procedures yielded confidence in the method proposed as they aligned well with findings from other studies, while at the same time, signified challenges inherent in multi-year crop phenological development modeling (influence of natural and anthropogenic factors).

Models Development and Evaluation
The main source of single year model misclassifications was due to similarities in the NDVI-based phenology curves. These similarities in the NDVI curves, representing different crops, can be potentially attributed to parallel farming operation schedules coupled with comparable crop biomass changes over the growing cycle. Challenges in signal separability between the NDVI-based phenology curves of major row crops have also been reported in other studies [26,32,47]. These similarities in temporal signatures pose a challenge to classifier algorithms, although in this study, it was more pronounced in the 2000 than in the 2005 dataset. dimensionality standardization (piecewise cubic polynomial). However, these datasets are not immune to noise. Inter-field variations could be potentially attributed to gaps in the dataset (clouds) incorrectly filled by subsequent steps, inconsistencies in the CDL layer used as the reference, or contaminated pixels (cloud shadows) affecting the individual time period. This uncertainty is often addressed with increased temporal resolution, but remains a challenge when the objective is to classify historic datasets where temporal resolution is limited. Annual variations also influenced the model's robustness, represented by the ability to generalize when presented with datasets not used during the training stage. Single year developed models were vulnerable to anthropogenic and natural induced watershed-wide changes. Anthropogenic factors such as crop conversion influenced how well a particular crop type class is represented during the model development (training stage). Based on reference data generated from CDLs from 2000 to 2005, an increase in corn/maize (from 182 to 220 fields) and soybean (from 292 to 676 fields) was observed in this watershed, but a decrease in spring wheat (from 534 to 377 fields) and sunflower (from 47 to 12 fields). This regional trend is consistent with national crop conversion moves toward corn/maize and soybean rotations [49,50]. Natural factors are simply due to climatic variations [48]. Developing models based on multiyear datasets has addressed uncertainties by annual variations. Results from the 2000-2005 model are in agreement with similar studies in China [27], Brazil [31], and Kansas, USA [26], despite the fact these studies used a higher temporal resolution either by using MODIS or combining different platforms.
It is important to acknowledge the accuracy of CDL layers calculated at the pixel scale.

Model Application
The lack of historic reference data prevented the absolute accuracy assessment of crop type generated for 1985, 1990, and 1995. However, relative comparisons between the different model's predictions demonstrated agreement in the dominant crop area estimation and between pairs of models. The agreement between models can also be observed in the crop-type conversion trend (Figure 8). Spatial models ( Figure 6) in conjunction with temporal models (Figure 8) support the development of input databases capturing the dominant crop-conversion trends, which is important information for long-term hydrological/climatic studies designed to characterize existing conditions over time. Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 16

Conclusions
Long-term hydrological and/or climatological modeling studies of agriculture watersheds depend on reliable crop-specific land cover information. This study proposed and evaluated methods for deriving such datasets solely based on the time series of satellite imagery through integrating data processing methods (filtering and data gap filling) with non-parametric classifying algorithms (random forest algorithm). Uncertainties in mis-registration between images was addressed by working in the object domain, represented by farm fields, rather than the pixel domain. Additional advantages of working in the object domain versus the pixel domain include the possibility of increasing the temporal resolution by being able to use objects that are partially contaminated (either clouds or no-data pixels), the significant reduction in data volume (from millions to thousands), and the generation of smoothed datasets (field averaged). In contrast, field boundaries were assumed to be static over time. In a small number of fields, the cultivation of more than one crop type in a single season (field is split into two or three crops) was observed. This was addressed by developing custom zonal statistical algorithms to consider only fields where a single crop-type covered more than 50% of the field. Future studies could investigate the use of image segmentation algorithms as possible alternatives to develop dynamic image objects (field boundaries changing from year to year).
Predictive model development based on single year datasets (NDVI and CDL) yielded adequate overall accuracy results for both the training and testing datasets. The main sources of uncertainties in the model development were phenological similarities between the crop-type and representativeness of crop type within the training datasets. The later uncertainty was addressed through the utilization of the multi-year training dataset. This approach has the advantage of increasing crop type representation during the training phase and capturing potential yearly variations due to natural and anthropogenic factors.
The integration of machine learning algorithms, enhanced signal filtering, and dimensionality standardization through function fitting with remotely sensed time series of spectral indices to develop historic crop data layers at the field scale is a feasible alternative. This approach has the potential to temporally extend CDL coverage in the USA, and therefore support long-term hydrological/climatological modeling of agriculture watersheds. Accurate crop-type characterization influences how farming practices are modeled in time and space and therefore how non-point source pollutants are produced, transported, and deposited; thus, supporting the development of effective mitigation strategies.

Conclusions
Long-term hydrological and/or climatological modeling studies of agriculture watersheds depend on reliable crop-specific land cover information. This study proposed and evaluated methods for deriving such datasets solely based on the time series of satellite imagery through integrating data processing methods (filtering and data gap filling) with non-parametric classifying algorithms (random forest algorithm). Uncertainties in mis-registration between images was addressed by working in the object domain, represented by farm fields, rather than the pixel domain. Additional advantages of working in the object domain versus the pixel domain include the possibility of increasing the temporal resolution by being able to use objects that are partially contaminated (either clouds or no-data pixels), the significant reduction in data volume (from millions to thousands), and the generation of smoothed datasets (field averaged). In contrast, field boundaries were assumed to be static over time. In a small number of fields, the cultivation of more than one crop type in a single season (field is split into two or three crops) was observed. This was addressed by developing custom zonal statistical algorithms to consider only fields where a single crop-type covered more than 50% of the field. Future studies could investigate the use of image segmentation algorithms as possible alternatives to develop dynamic image objects (field boundaries changing from year to year).
Predictive model development based on single year datasets (NDVI and CDL) yielded adequate overall accuracy results for both the training and testing datasets. The main sources of uncertainties in the model development were phenological similarities between the crop-type and representativeness of crop type within the training datasets. The later uncertainty was addressed through the utilization of the multi-year training dataset. This approach has the advantage of increasing crop type representation during the training phase and capturing potential yearly variations due to natural and anthropogenic factors.
The integration of machine learning algorithms, enhanced signal filtering, and dimensionality standardization through function fitting with remotely sensed time series of spectral indices to develop historic crop data layers at the field scale is a feasible alternative. This approach has the potential to temporally extend CDL coverage in the USA, and therefore support long-term hydrological/climatological modeling of agriculture watersheds. Accurate crop-type characterization influences how farming practices are modeled in time and space and therefore how non-point source pollutants are produced, transported, and deposited; thus, supporting the development of effective mitigation strategies.