Fusion of Moderate Resolution Earth Observations for Operational Crop Type Mapping

Crop type inventory and within season estimates at moderate (<30 m) resolution have been elusive in many regions due to the lack of temporal frequency, clouds, and restrictive data policies. New opportunities exist from the operational fusion of Landsat-8 Operational Land Imager (OLI), Sentinel-2 (A & B), and Sentinel-1 (A & B) which provide more frequent open access observations now that these satellites are fully operating. The overarching goal of this research application was to compare Harmonized Landsat-8 Sentinel-2 (HLS), Sentinel-1 (S1), and combined radar and optical data in an operational, near-real-time (within 24 h) context. We evaluated the ability of these Earth observations (EO) across major crops in four case study regions in United States (US) production hot spots. Hindcast time series combinations of these EO were fed into random forest classifiers trained with crop cover type information from the Cropland Data Layer (CDL) and ancillary ground truth. The outcomes show HLS achieved high (>85%) accuracies and the ability to provide insight on crop location and extent within the crop season. HLS fused with S1 had, at times, a higher accuracy (5–10% relative overall accuracy and kappa increases) within season although the combination of fused data was minimal at times, crop dependent, and the accuracies tended to converge by harvest. In cloud prone regions and certain temporal periods, S1 performed well overall. The growth in the availability of time dense moderate resolution data streams and different sensitivities of optical and radar data provide a mechanism for within season crop mapping and area estimates that can help improve food security.


Introduction
Assessment and monitoring of crop type and extent is one of the most critical information needs for food security followed by indicators of health and production.Stakeholders at all levels require reliable and timely information to make sound decisions that can improve inventory, investments, and mitigations.Since the 1970's, projects such as the Large Area Crop Inventory Experiment (LACIE) and Agricultural and Resources Inventory Surveys Through Aerospace Remote Sensing (AgRISTARS) have been leveraging spaceborne Earth observations to map crop type and extent.These early works have evolved into many national scale systems and decision support tools that now leverage moderate resolution optical data for delineating croplands.For example, in the United States (US) the Department of Agriculture, National Agricultural Statistics Service (USDA-NASS) generates annually the Cropland Data Layer (CDL).In the past five years, CDL inputs have focused on using Landsat along with inputs from Deimos-1, UK-Disaster Monitoring Constellation (DMC), and Sentinel-2 to generate nominal 30 m pixel resolution maps of crop type.Within the season they are used internally by NASS to supplement planted area survey information but after the season they are released to the public and as such have been used in a variety of landcover applications.
The increase in moderate resolution sensors, open and free access policies, and operational and systematic coverage have created a new era of opportunities to improve operational crop type assessment and monitoring (that is, Reference [1]).In particular, one promising source for driving next-generation products is the potential of near-daily, moderate resolution, multispectral optical data with the harmonization of Sentinel-2A and B with Landsat-8 Operational Land Imager (OLI).For the first time, the applied science community has access to systematic, near-daily, moderate resolution optical data.These scales of multispectral time series can potentially enable the transferability of approaches that have been developed using sensors such as Moderate Resolution Imaging Spectroradiometer (MODIS) (for example, References [2][3][4][5]) to support food security.Recently, the globes first open access 30 m resolution crop area map that leveraged Landsat and cloud computing was released as part of the Global Food Security-Support Analysis Data @ 30-m (GFSAD30) Project [6].Other example efforts are showing the power of blending these optical sensors together, for example, mapping national scale within season estimates of soybeans in the US [7].
Concordantly, the growth in the availability of Synthetic Aperture Radar (SAR) from Sentinel-1A and B has created new opportunities for operational crop type mapping since the launch in 2014.The all-weather capability, active sensing system that operates independent of cloud cover and sun illumination, and sensitivity to surface and subsurface characteristics (that is, dielectric constant, roughness, orientation) make SAR particularly useful for mapping agricultural and field characteristics.However, historically SAR applications for crop monitoring have been much less relative to optical data such as Landsat or MODIS.The reasons for this include, limited availability and little to no open access operational observations, no consistent, large-area acquisition strategies, poor quality digital elevation models required for processing, and complex data structures relative to optical data.
A few institutions, such as Agriculture and Agri-Food Canada, have been using SAR for operational crop inventory, benefiting from access to ample Radarsat [8].However, while many SAR sensors (that is, ERS-1, ENVISAT ASAR, TerraSAR-X, Radarsat, ALOS) have been utilized for research on crop mapping [9][10][11][12][13], few programs use SAR for operational crop mapping.Rice crops, given their role in global food security and the tendency for cultivation in cloud prone regions, have long utilized SAR for monitoring.The Asian Rice Crop Estimation and Monitoring (Asia-RiCE) and Group on Earth Observations Global Agricultural Monitoring (GEOGLAM) have spearheaded the use of rice mapping using SAR and are beginning to scale up large-area operational systems for rice production monitoring [14,15].Over the next few years, the Joint Experiment for Crop Assessment and Monitoring (JECAM) will further lead an initiative across global sites to develop SAR-optical crop type mapping techniques as part of the GEO Agricultural Monitoring Community of Practice.
The overarching goal of this research application was to evaluate operational fusion of Sentinel-1 with Sentinel-2 and Landsat-8 for near real-time crop type mapping.Only very recently has the science community had the opportunity to pragmatically blend these moderate resolution optical and SAR data for crop inventory in an operational context.To achieve operational fusion of these moderate resolution sensors for crop type mapping, the strengths and limitations need to be identified as well as the development of methodologies.Objectives include the use of dense time series moderate resolution EO, determination of the accuracy of sensors individually and combined, the performance for a given crop, and when acceptable accuracy levels occur within the crop season.

Study Areas
Four study regions approximately 110 × 110 km each were selected as case study sites (Figure 1).These regions were selected to represent major agricultural production hot spots in the USA with active field research that cover diverse bioclimates, soils, calendars, cloud probabilities, and management practices to provide a robust set of test landscapes.These included northwest Ohio (NWO), northeast Arkansas (NEA), southeast South Dakota (SED), and the northern Sacramento Valley, California (VAI).Major crops (>5% of landscape) in NWO and SED include corn, soybeans, winter wheat, and pasturelands.NEA includes cotton and rice production in addition to corn and soy.VAI has a diverse crop matrix with areas of grapes, tomatoes, corn, rice, alfalfa, sunflower, clover, almonds, and walnuts among other specialty crops (for example, onions, peas, watermelon, carrots).Rain season, crop calendars, and management practices vary within and among the study regions creating a robust set of calibration and validation crop type sites.

Study Areas
Four study regions approximately 110 × 110 km each were selected as case study sites (Figure 1).These regions were selected to represent major agricultural production hot spots in the USA with active field research that cover diverse bioclimates, soils, calendars, cloud probabilities, and management practices to provide a robust set of test landscapes.These included northwest Ohio (NWO), northeast Arkansas (NEA), southeast South Dakota (SED), and the northern Sacramento Valley, California (VAI).Major crops (>5% of landscape) in NWO and SED include corn, soybeans, winter wheat, and pasturelands.NEA includes cotton and rice production in addition to corn and soy.VAI has a diverse crop matrix with areas of grapes, tomatoes, corn, rice, alfalfa, sunflower, clover, almonds, and walnuts among other specialty crops (for example, onions, peas, watermelon, carrots).Rain season, crop calendars, and management practices vary within and among the study regions creating a robust set of calibration and validation crop type sites.

Sentinel-1 Processing
Dual polarization (VV + VH) Sentinel-1 Interferometric Wide (IW) Ground Range Detection (GRD) data were obtained from the Alaska Satellite Facility (ASF) for this application.Sentinel-1A carries a C-band imager at 5.405 GHz with an incidence angle between 20° and 45°.The platform follows a Sun-synchronous, near-polar, circular orbit at a height of 693 km.These C-band radar data were first transformed into Beta Naught to assist in radiometric terrain correction and thermal noise removal using the look-up table (LUT) within the metadata.Radiometric distortion caused by the terrain slope was corrected using the method proposed by Reference [16] to achieve Gamma Naught γ°, which is less dependent of the incidence angle than the Sigma Naught σ°.The inherent speckle noise was filtered using a Boxcar filter method with a 5-by-5 window size.Finally, with the assistance of the precise orbit and 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM), all data were geocoded using the range-Doppler approach with the output resolution of 30 by 30 m (to match a nominal application scale in this case).All radar preprocessing were executed using custom open access routines we built using Python, C, and GDAL and are available for sharing with author contact.At the time of this application, only Sentinel-1A was uniformly operational over these study domains.

Sentinel-1 Processing
Dual polarization (VV + VH) Sentinel-1 Interferometric Wide (IW) Ground Range Detection (GRD) data were obtained from the Alaska Satellite Facility (ASF) for this application.Sentinel-1A carries a C-band imager at 5.405 GHz with an incidence angle between 20 • and 45 • .The platform follows a Sun-synchronous, near-polar, circular orbit at a height of 693 km.These C-band radar data were first transformed into Beta Naught to assist in radiometric terrain correction and thermal noise removal using the look-up table (LUT) within the metadata.Radiometric distortion caused by the terrain slope was corrected using the method proposed by Reference [16] to achieve Gamma Naught γ • , which is less dependent of the incidence angle than the Sigma Naught σ • .The inherent speckle noise was filtered using a Boxcar filter method with a 5-by-5 window size.Finally, with the assistance of the precise orbit and 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM), all data were geocoded using the range-Doppler approach with the output resolution of 30 by 30 m (to match a nominal application scale in this case).All radar preprocessing were executed using custom open access routines we built using Python, C, and GDAL and are available for sharing with author contact.At the time of this application, only Sentinel-1A was uniformly operational over these study domains.

Harmonized Landsat-8 Sentinel-2
The Harmonized Landsat-8 Sentinel-2 (HLS) is an ongoing processing chain effort to take advantage of widely accessible moderate resolution optical platforms to increase the potential number of temporal revisits (hls.gsfc.nasa.gov).This processing chain is evaluated over validation sites from AERosol RObotic NETwork, Fluxnet, Societal Applications in Fisheries & Aquaculture using Remotely Sensed Imagery, and Baseline Surface Radiation Network programs.The harmonization approach (1) grids imagery to a common pixel resolution, projection, and spatial extent (that is, tile); (2) applies atmospheric correction, the Landsat 8 Surface Reflectance Code (LaSRC) processor [17,18], relying of on a long heritage of Landsat processing; (3) applies LaSRC cloud masking using efficient cloud and cloud-shadow modules for Landsat-8 data [17,18] while Fmask [19] is used for Sentinel-2; (4) adjusts to represent the response from a common spectral bandpass and BRDF-normalization keying off a single, global and constant BRDF shape that produces satisfying BRDF normalization over a limited range of view zenith angle [20,21]; (5) perform geographic co-registration using the Automated Registration and Orthorectification Package (AROP) [22] to warp and co-register data.Note that the bandpass adjustment corresponds to a linear fit between the equivalent spectral bands.The regression coefficients were computed based on 500 hyperspectral spectra selected on 160 Hyperion scenes globally distributed.At the end of this harmonization processing chain, a temporal frequency of 1-3 days using synchronized platforms is achieved (hls.gsfc.nasa.gov).In this application reflectances, common indices (that is, NDVI, LSWI) that have shown value for classifying crops (that is, [23][24][25][26]) and thermal bands (TIR) were included as inputs into the classifier.

Classification and Assessment
Summarizing, S1, HLS, and blended HLS + S1 stacks were used in a hindcast approach (Figure 2, Table 1) to classify major crop types in the four case study regions.For a classifier, we used the ensemble, machine-learning, Random Forest (RF) algorithm [27].Random forest is a flexible and powerful nonparametric technique that many mapping applications have recently implemented for a range of studies (for example, [28][29][30][31]).A random forest is generated through the creation of a series of classification/decision trees using bootstrapping, or resampling with replacement.Tuning parameters, such as the number of trees and the number of split candidate predictors, are generally chosen based on the out-of-bag (OOB) prediction error.In this application, we implement the random forest classifier from python sklearn and tuned with the exception of the number of trees.The number of trees is set to a large number (200) because compared with other classifiers such as Support Vector Machine (SVM), the RF uses out-of-bag (OOB) samples for cross-validation, and once the OOB errors stabilize at a sufficiently large number of trees, the training can be concluded.The number of m variables that are selected randomly as candidates for splitting is set to the default which is the square root of the number of total input predictors p.It should be noted that when the number of variables p is large, but the fraction of relevant variables is small, RF can perform poorly with a small m.We further use a variable importance metric based on the Gini index, where predictors with higher importance are used more often to create a split, to understand how different days, sensors and bands relate to model performance, noting that varying scales of measurement between the predictors may also influence the importance ranking [32,33].
Due to the diverse imaging modes, acquisition dates, and the invalid values caused by cloud within HLS, we interpolated the images stacks to standardized 15-day time windows between the day of year 96 (6 April) and 308 (4 November).A 15-day interval was chosen to ensure at least one observation from each of the sensors was available in every time window and further supported by prior knowledge and subjective interpretation of plots (Figure 3).The interpolation was carried out using a multi-temporal time series routine that allows one to explore the temporal domain of a large space-time data set in an efficient manner.The routine, by reading small spatial chunks of all available time layers, carries out spatial location specific calculations on all-time layers grouped by time (that is, the year) and outputs maps of the time calculated values.Calculation algorithms are contained in separate modules that can be strung together in a user-specified step order for more complex calculations.In this application, we use a two-module step algorithm.It begins with a recomposite step, which calculates a band specific average for a 15-day window at each pixel and concludes with a gapfill step, which linearly interpolates the 15-day band specific average to fill in periods with no observations.complex calculations.In this application, we use a two-module step algorithm.It begins with a recomposite step, which calculates a band specific average for a 15-day window at each pixel and concludes with a gapfill step, which linearly interpolates the 15-day band specific average to fill in periods with no observations.We qualitatively considered various options for smoothing (that is, References [34,35]) following the best practices and applications [36][37][38][39].Ultimately, the focus of this work was blending the optical and radar data for crop type classification case study; thus, a linear interpolation was sufficient and simple to execute.If the first 15-day window was missing for a specific pixel and band, the spatial average of that band for all crop areas was used to fill the missing value.The multi-temporal processing for one single example pixel is shown in Figure 3.The advantage of this processing is an extremely efficient implementation and its flexibility to be used in a near real-time application since the approach gapfill interpolates using only past dates.Further investigation is needed to assess the uncertainties added by this approach, in particular, when few time points are available overall.
Crop type training and validation came from the 2016 USDA NASS CDL.Details on the creation of CDL can be found in Reference [40].The CDL is a very good proxy for ground information because it provides pixel-level accuracies of around 95% for major commodities in intensively farmed regions (Table 2).For training, crop stratified random pixel (30 m resolution) samples of size 1000 in each of the four case study areas were selected from a sampling unit.Samples were restricted to homogenous patches of 5 by 5 windows with the same CDL class and pixels having at least five Landsat 8, three Sentinel 2 and five Sentinel 1 valid overpass dates during the crop season (DOY 96 to 308).The subsample locations were used to extract time series' from multitemporal processed raster stacks of Harmonized Landsat and Sentinel 2 bands as well as Sentinel 1 bands VV and VH.Validation was carried out using all CDL major crop type pixels over each study area.Error matrices were generated based on these outcomes.We qualitatively considered various options for smoothing (that is, References [34,35]) following the best practices and applications [36][37][38][39].Ultimately, the focus of this work was blending the optical and radar data for crop type classification case study; thus, a linear interpolation was sufficient and simple to execute.If the first 15-day window was missing for a specific pixel and band, the spatial average of that band for all crop areas was used to fill the missing value.The multi-temporal processing for one single example pixel is shown in Figure 3.The advantage of this processing is an extremely efficient implementation and its flexibility to be used in a near real-time application since the approach gapfill interpolates using only past dates.Further investigation is needed to assess the uncertainties added by this approach, in particular, when few time points are available overall.
Crop type training and validation came from the 2016 USDA NASS CDL.Details on the creation of CDL can be found in Reference [40].The CDL is a very good proxy for ground information because it provides pixel-level accuracies of around 95% for major commodities in intensively farmed regions (Table 2).For training, crop stratified random pixel (30 m resolution) samples of size 1000 in each of the four case study areas were selected from a sampling unit.Samples were restricted to homogenous patches of 5 by 5 windows with the same CDL class and pixels having at least five Landsat 8, three Sentinel 2 and five Sentinel 1 valid overpass dates during the crop season (DOY 96 to 308).The subsample locations were used to extract time series' from multitemporal processed raster stacks of Harmonized Landsat and Sentinel 2 bands as well as Sentinel 1 bands VV and VH.Validation was carried out using all CDL major crop type pixels over each study area.Error matrices were generated based on these outcomes.We further qualitatively assessed training data using a limited set of geofield photos taken during the crop season.These field-level "windshield drive-by" photos were collected using a smartphone app.All the geofield photos are linked to shape files or keyhole markup language (KML) files to store, display, and share photos.These photos are available for viewing and sharing at www.eomf.ou.edu/ photos.Further, we used a limited set of Common Land Unit (CLU) polygons from the USDA Farm Service Agency to generate mean summary statistics from the integrated data stacks.These CLU data were also qualitatively used to compare training data statistically.Together, these ultimately provide high confidence in the use of the calibration and validation data.

Results
The cloud masks applied to the HLS imagery removed an average of 35, 37, 59, and 29% of valid pixels for NEA, NWO, SED, and VAI, respectively, across time (DOY 90-310) (Figure 4).Noted is that we did not perform regional tuning of the cloud mask algorithms given the operational context of this research application.The heat maps illustrating the number of quality observations shows that the valid counts range from 4-40 across the four case study sites.Further, the number of quality observations varies within the study regions due to the path row footprints, cloud, ascending versus descending orbits and the effects of the observations strategies.The highest frequency of observations is apparent in the path overlap regions.
Accuracy measures range by region, sensor, time, and crop type.Note these error matrices reflect the model classes in that only major agricultural land uses were used (that is, no built, water, forested, and so forth), which, therefore, influences the accuracy and kappa outcomes.Using NEA (Figure 5) as an example, S1-only overall accuracies (OA) and kappa values for the first (DOY 96) period were 30% and 7, respectively, while the outcomes at the end of the season (DOY 300) were 82% and 72, respectively (Table 3).For HLS, the overall accuracies and kappa values for the first period were 47% and 27 while the outcomes at the end of the season were 92% and 87, respectively.For integrated stacks (HLS + S1) the overall and kappa outcomes were 48% and 28 initially while the end of the season values were 94% and 90, respectively.Users and Producers Accuracies (UA/PA) for the first and end of season periods for NEA for major crops including corn, cotton, rice, and soybean are illustrated in Table 3. Corn had the most confusion with S1-only in NEA with UA of 52% by end of the season compared to 80% when using HLS.In VAI, almonds had an integrated UA of 27 and fallow cropland and winter wheat had UA slightly below 60% while all other crops had UA and PA over 60%, including corn, rice, alfalfa, tomatoes, grapes, grassland, and walnuts (Figure 6A,B).To succinctly summarize the hundreds of classifications (sensors × time × regions), Figure 7 summarizes the overall accuracy and kappa values for each site over time and by sensor type.Evident is the increasing accuracies as seasonal crops develop.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17 Users and Producers Accuracies (UA/PA) for the first and end of season periods for NEA for major crops including corn, cotton, rice, and soybean are illustrated in Table 3. Corn had the most confusion with S1-only in NEA with UA of 52% by end of the season compared to 80% when using HLS.In VAI, almonds had an integrated UA of 27 and fallow cropland and winter wheat had UA slightly below 60% while all other crops had UA and PA over 60%, including corn, rice, alfalfa, tomatoes, grapes, grassland, and walnuts (Figure 6A,B).To succinctly summarize the hundreds of classifications (sensors × time × regions), Figure 7 summarizes the overall accuracy and kappa values for each site over time and by sensor type.Evident is the increasing accuracies as seasonal crops develop.The variable importance was measured using the Gini index.Looking at the single highest Gini index by model, a given region typically maxes out using a certain variable (band) on a certain day of the year.These outcomes coincide with the overall accuracy and kappa plateaus for a given region.For NWO and SED, this is 197, while for NEA, this is 239 and in VAI it depends on the sensor (that is, HLS NDVI 113) or data availability.In VAI, there is a visible influence of data by sensor availability, which also factors in which bands have the highest Gini index.In VAI, NDVI, LSWI, and S2, the red edge-1 are intermixed over time as the highest ranked.In SED, NDVI followed by S2 red edge-1 show up as the most important from DOY 141 to 239.In NWO and NEA, LSWI, VV, and SWIR are ranked the highest.Cirrus and TIR were often ranked the least important across the four case study sites and over time.

Discussion
Generally, the HLS performed very well in these case studies of mapping major crop type in four US production hot spots.VAI with its 10 diverse crops was the only site which did not achieve >85% overall accuracy keeping in mind outcomes are relative (only major crop classes considered) and should be taken in the context of the study design.The S1 alone performed adequately while the combination of the radar and optical satellite imagery tended to out-perform any individual sensor  The variable importance was measured using the Gini index.Looking at the single highest Gini index by model, a given region typically maxes out using a certain variable (band) on a certain day of the year.These outcomes coincide with the overall accuracy and kappa plateaus for a given region.For NWO and SED, this is 197, while for NEA, this is 239 and in VAI it depends on the sensor (that is, HLS NDVI 113) or data availability.In VAI, there is a visible influence of data by sensor availability, which also factors in which bands have the highest Gini index.In VAI, NDVI, LSWI, and S2, the red edge-1 are intermixed over time as the highest ranked.In SED, NDVI followed by S2 red edge-1 show up as the most important from DOY 141 to 239.In NWO and NEA, LSWI, VV, and SWIR are ranked the highest.Cirrus and TIR were often ranked the least important across the four case study sites and over time.

Discussion
Generally, the HLS performed very well in these case studies of mapping major crop type in four US production hot spots.VAI with its 10 diverse crops was the only site which did not achieve >85% overall accuracy keeping in mind outcomes are relative (only major crop classes considered) and should be taken in the context of the study design.The S1 alone performed adequately while the combination of the radar and optical satellite imagery tended to out-perform any individual sensor earlier in the growing season.However, in many cases, the fusion of optical and radar only provide incremental improvements.For any given region, major crop, or time of year example arguments can be made of which sensor or combination is superior.We emphasize that this is a relatively simple case study testing ability to distinguish major crop types (categorical variable) only with the ensemble machine learning classifier.
We note that CDL was used as the primary calibration and validation data in this research application and this has inherent uncertainties.With the robust sampling approach described in combination with limited geofield photos and CLU considerations, we feel confident in the quality of the calibration and validation.The CDL PA and UA all have crop specific accuracies of over 93%, excluding corn, winter wheat, and pistachios in VAI.Therefore, it is possible these uncertainties influenced the outcomes of VAI.However, in our experience, CLUs, which are used to train CDL, have embedded errors and the misclassifications in CDL often tend to be related to noise (edge effect, speckling) so we suggest that this is as robust an overall approach as operationally feasible.
Clouds are often cited as the main advantage in using radar.Here the cloud masks generated with the HLS products were not further tuned and taken and applied as is.In the Southeast South Dakota study region, 38 HLS overpasses had more than 90% of the pixels flagged.It is very likely that these numbers could be lowed (improved) with regional tuning and cloud mask refinement.Masking clouds remains an evolving and active research topic [41,42].However, the classification outcomes indicate the multitemporal routine, of compositing and interpolating the quality HLS observations, and the classification approach still out-performed the S1 only data by 5% and 8 for overall accuracy and kappa, respectively.Potentially, HLS in geographies like South or Southeast Asia with persistent clouds might be severely limited during crucial periods for food security decision making.These outcomes are also shaped by the study design and, in fact, only two bands are available from S1 IW mode (that is, VV, VH) versus several bands available from HLS and the use of random forest.It is possible that techniques leveraging quad polarization or repeat-pass interferometric information with more sensitivity tuned to the number and density, dielectric contrast, orientation and shape, or surface roughness of targets would improve radar only accuracy (for example, References [43,44]).
The timing of quality observations had an influence on the outcomes.The OA and kappa tend to plateau around DOY 200 for southeast Dakota and northwest Ohio with maximum accuracy leveling out.This is likely due to the relatively simple crop classification schemes of only corn, soy, and grasses (winter wheat and/or pasture).Northeast Arkansas did not plateau until DOY 240, and generally, S1 tended to plateau a few 15-day periods later than HLS, with the classifier using the additional observations to statistically differentiate between categories.This further amplifies the value of the additional parameters (bands) in terms of random forest computational needs.On the other hand, this approach allowed for a near real-time implementation.
The timing and number observations are pronounced in the variable importance analysis.The 'most valuable' variable, according to this study context, can shift throughout the time series and does not necessarily have biogeophysical meaning tied to crop phenological attributes (for example, chlorophyll, biomass, leaf area index, roughness).Since we fit models iteratively in time there are complex interactions between the day of year and sensor band.In particular, the band with the highest Gini index for a particular image day input changes depending on the other dates available in the model.In Figure 8, we show this for each of the by-day NEA models when all HLS and S1A bands are included.For example, the band with the highest Gini index for the first image recomposite period (DOY 90) includes S2 Rededge-1, TIRS1, NIR, and NDVI.However, overall, the VH, LSWI, NIR, and SWIR bands populated the top variables during the key early season periods of crop development (Figure 8).This likely indicates that volume, leaf moisture or equivalent water thickness, green biomass vigor, and residue and organic matter are driving these parameters to the top of the classifier importance rankings.Parameters (bands) that showed up as the least important were cirrus and TIRS with the earlier periods (for example, DOY 90) getting the least weight.
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 17 that these numbers could be lowed (improved) with regional tuning and cloud mask refinement.Masking clouds remains an evolving and active research topic [41,42].However, the classification outcomes indicate the multitemporal routine, of compositing and interpolating the quality HLS observations, and the classification approach still out-performed the S1 only data by 5% and 8 for overall accuracy and kappa, respectively.Potentially, HLS in geographies like South or Southeast Asia with persistent clouds might be severely limited during crucial periods for food security decision making.These outcomes are also shaped by the study design and, in fact, only two bands are available from S1 IW mode (that is, VV, VH) versus several bands available from HLS and the use of random forest.It is possible that techniques leveraging quad polarization or repeat-pass interferometric information with more sensitivity tuned to the number and density, dielectric contrast, orientation and shape, or surface roughness of targets would improve radar only accuracy (for example, References [43,44]).The timing of quality observations had an influence on the outcomes.The OA and kappa tend to plateau around DOY 200 for southeast Dakota and northwest Ohio with maximum accuracy leveling out.This is likely due to the relatively simple crop classification schemes of only corn, soy, and grasses (winter wheat and/or pasture).Northeast Arkansas did not plateau until DOY 240, and generally, S1 tended to plateau a few 15-day periods later than HLS, with the classifier using the additional observations to statistically differentiate between categories.This further amplifies the

Conclusions
Overall, the outcomes of these case study crop type classifications highlight the utility of higher temporal frequency moderate resolution data streams.Both HLS and S1 offer contributions to characterizing major crop types with the combination of radar and optical adding value (5-10% relative accuracy increases earlier and during the crop season).By the end of the crop season, the overall accuracy and kappa tend to converge showing equivalent capability in these regions.The harmonization of Landsat-8 and Sentinel-2 performed best and evident is the ability to distinguish major crop types within a season with high accuracy approaching harvest.Typically, observations from these sensors are available within 12-24 h providing a mechanism for near-real-time capabilities.However, it should be noted that currently HLS processing is not on demand and requires additional harmonization processing.In regions with chronic and persistent cloud cover, this work indicates Sentinel-1 can be very effective in providing major crop type inventory estimates within the season.With ESA striving for 24-h latency on S1 this should be able to support operational production.Given the short latency and quality calibration of these sensors, the science community is now beginning to take advantage of these opportunities to drive decision support tools and inventory programs to enhance food security with within season estimates of extent.Being able to monitor and assess production at field scale will drastically enhance decision making.Future work can consider evaluating fusion in other regions and crops, transferability of models across space and time (years), using additional ground truth for accuracy and uncertainty quantification, and scaling over large regions operationally are potential next steps.As more moderate resolution sensors come online in an operational context the agricultural monitoring community can potentially support within season major crop type production estimates.This can help to meet the most important need for food security.

Figure 1 .
Figure 1.The four diverse case study sites located in major production hot spots in the USA with varying bioclimates, calendars, and managements.

Figure 1 .
Figure 1.The four diverse case study sites located in major production hot spots in the USA with varying bioclimates, calendars, and managements.

17 Figure 2 .
Figure 2. The conceptual workflow fusing multitemporal radar and optical moderate resolution imagery.

Figure 2 .
Figure 2. The conceptual workflow fusing multitemporal radar and optical moderate resolution imagery.

Figure 3 .
Figure 3.An example (A) time-frequency collected by sensor type across the northeast Arkansas (NEA) case study region.An example (B) multitemporal time series routine post processing Sentinel-1 VH observations a single corn pixel from NEA. Open circles (raw) give the VH backscatter γ° from the original image overpass date (2016).Blue lines (recomposite) are the mean value for each 15-day window shown by the vertical gray dotted lines, and the red line illustrates the final smoothed/interpolated classification inputs.

Figure 3 .
Figure 3.An example (A) time-frequency collected by sensor type across the northeast Arkansas (NEA) case study region.An example (B) multitemporal time series routine post processing Sentinel-1 VH observations a single corn pixel from NEA. Open circles (raw) give the VH backscatter γ • from the original image overpass date (2016).Blue lines (recomposite) are the mean value for each 15-day window shown by the vertical gray dotted lines, and the red line illustrates the final smoothed/interpolated classification inputs.

Figure 4 .
Figure 4.The heat maps by sensor and region showing cloud masked observation frequency.

Figure 5 .
Figure 5. Example end of season crop type classification maps in northeast Arkansas by sensors showing a generally strong ability to map major crop type within season at the field scale.The highlighted are example S1 that are only misclassification among corn, soy, and cotton.

Figure 4 .
Figure 4.The heat maps by sensor and region showing cloud masked observation frequency.

Figure 4 .
Figure 4.The heat maps by sensor and region showing cloud masked observation frequency.

Figure 5 .
Figure 5. Example end of season crop type classification maps in northeast Arkansas by sensors showing a generally strong ability to map major crop type within season at the field scale.The highlighted are example S1 that are only misclassification among corn, soy, and cotton.

Figure 5 .
Figure 5. Example end of season crop type classification maps in northeast Arkansas by sensors showing a generally strong ability to map major crop type within season at the field scale.The highlighted are example S1 that are only misclassification among corn, soy, and cotton.

Figure 6 .
Figure 6.(A) The end of season crop type products for northwest Ohio, southeast South Dakota, and Sacramento Valley California by sensor type and combination; zoom-in continued on (B).(B) Representative zoomed in snapshots from (A), the end of season crop type products for southeast South Dakota, northwest Ohio, and Sacramento Valley California by sensor type and combination.

BFigure 6 . 17 Figure 7 .
Figure 6.(A) The end of season crop type products for northwest Ohio, southeast South Dakota, and Sacramento Valley California by sensor type and combination; zoom-in continued on (B).(B) Representative zoomed in snapshots from (A), the end of season crop type products for southeast South Dakota, northwest Ohio, and Sacramento Valley California by sensor type and combination.Remote Sens. 2018, 10, x FOR PEER REVIEW

Figure 7 .
Figure 7.The overall accuracy and kappa values by sensor types over the crop season in the four case study regions.

Figure 8 .
Figure 8.The summary of the highest ranked parameter (band) across the four case study regions over time.For example, on DOY 127 following bands can be seen as highest ranked for the given regions: NEA VH and NIR-HLS; NWO SWIR1-HLS, NDVI-HLS, and VV; SED NDVI-HLS and LSWI-HLS; and VAI NIR-HLS, NDVI-HLS, and LSWI-HLS.

Figure 8 .
Figure 8.The summary of the highest ranked parameter (band) across the four case study regions over time.For example, on DOY 127 following bands can be seen as highest ranked for the given regions: NEA VH and NIR-HLS; NWO SWIR1-HLS, NDVI-HLS, and VV; SED NDVI-HLS and LSWI-HLS; and VAI NIR-HLS, NDVI-HLS, and LSWI-HLS.

Table 1 .
The moderate resolution Earth observation inputs for crop classifier.Sensors Band Code Band Name Wavelength (Micrometers) CA Coastal Aerosol 0.44 B Blue 0.48

Table 1 .
The moderate resolution Earth observation inputs for crop classifier.

Table 2 .
The cropland data layer accuracies statistics for major crops across four case study regions.

Table 3 .
An example error matrix for one classification scheme in northeast Arkansas (NEA) by sensor combinations at Day Of Year 96 and 300 highlighting shifts in accuracy.