Knowledge of crop-specific planted area, conditions, and expected yields is critical for informing agricultural markets and decision making as well as ensuring food security globally. Official reports by agricultural ministries, such as the USDAâĂŹs Crop Progress and Condition Reports, are widely followed by market analysts and traders to inform decision-making around agricultural markets and futures. Thus, uncertainty in public information including official reports for major food-producing countries results in uncertainty and volatility in agricultural markets and prices [22] . This uncertainty can arise when they are shocks to the system, e.g., record flooding or the COVID-19 pandemic in the 2019 and 2020 growing seasons respectively. Market tensions and perceptions of increased uncertainty about agricultural production caused by policy decisions and reduced access to information influence market prices and decisions that can disrupt market stability and food security.
Earth observations (EO) data provide a valuable complementary source of information about global crop progress and conditions, enabling timely analysis at field scales throughout the growing season to supplement traditional data collection and reporting efforts. Crop type masks that provide classifications for crop-specific land cover are critical inputs for agricultural monitoring quantities including estimates of acreage, production, condition, and yield. Prior work on crop type classification has achieved promising results using time series information to differentiate spectrally-similar crop types. However, models trained to distinguish crop types based on temporal patterns may not generalize well to inputs that have shifted in the time domain, for example due to delayed growth or planting as a result of extreme weather events or climate change.
An example of this is the 2019 growing season in the midwestern United States. In spring 2019, heavy snowfall and rains caused unprecedented flooding across the midwestern states and forced farmers to plant crops 1-2 months later than normal [15] . While such extreme weather events are historically rare, they are becoming more frequent and exacerbate uncertainty and confusion for markets and decision-makers. Thus, for crop type classification models to be used operationally, it is important that they 1) be resilient to shifting growing seasons, and 2) provide predictions within the growing season when uncertainty around planting and production is highest.
In this paper, we present a method for crop type classification that uses time series inputs at key growth stages detected at the pixel level to mitigate the effects of domain shift caused by shifts in planting timelines. Additionally, this method enables in-season predictions of crop type at each detected growth stage. We demonstrated the performance of this method for mapping corn and soybean crops in a study region of the midwestern US in northern Illinois, and showed that our growth stage-normalized approach gives better classification accuracy than using fixed time series inputs. While this study is limited to crop type classification, our approach could be used for distinguishing other land cover types that exhibit unique phenological signatures, such as other plant or tree species, where time series classification approaches may be affected by seasonal shifts or climate change.
Most prior work on crop type classification incorporates spectral and temporal information in the input to a machine learning classifier. Random forests and decision trees are widely used for crop type classification due to their good performance in many case studies and their interpretability (e.g., [11, 25, 28, 31-33, 35-38, 40, 42, 43] ), though the input features extracted from the satellite observations vary. For example, Song et al. [35] used bagged decision trees to classify time series statistics extracted from optical Landsat spectral bands, whereas Wang et al. [43] used the coefficients of harmonic regression as input features for a random forest classifier. Kernel methods including support vector machines (SVMs) and kernelbased extreme learning machines [27] have also shown good performance for crop type classification (e.g., [9, 25, 31, 36, 37, 41, 42] )-for example, Feng et al. [9] extracted more than 100 features including spectral, textural, and phenological properties of optical observations for SVM classification.
One limitation of traditional machine learning methods such as random forests and SVMs is that they require input features to be extracted from the original data to achieve good performance, which may be sub-optimal for discriminating between target classes and often requires domain expertise. To address these limitations, recent work (e.g., [4-6, 10, 16, 19, 30, 45] ) has demonstrated improved classification performance using neural network-based approaches which automatically learn useful features for discriminating between target classes using a standard optimization procedure [21] . Zhong et al. [45] showed that a convolutional neural network (CNN) with 1D convolutions outperformed a long short term memory (LSTM) network [13] for classifying time series inputs of only one spectral index (enhanced vegetation index, or EVI). Cai et al. [5] showed that a deep neural network had better performance when information from multiple spectral bands and indices were included in the time series inputs. Kussul et al. [19] demonstrated better performance using an ensemble of 2D CNNs compared to 1D CNNs using time series from multiple optical bands and synthetic aperture radar (SAR) observations. Rußwurm et al. [30] showed that an LSTM gave better performance than a 2D-CNN for multispectral time series inputs that also included a neighborhood of pixels around the predicted pixel. Similarly, Garnot et al. [10] found that temporal features extracted in recurrent layers were more important than spatial features extracted in 2D convolutional layers for parcel-level crop type classification, but that hybrid methods combining recurrent and convolutional layers had the best performance. Brandt [4] showed that a network with capsule and attention layers in addition to CNN and LSTM layers performed better for crop type classification than a 2D CNN, LSTM, and combined CNN-LSTM. Following the success of 3D CNNs for modeling spatio-temporal relationships in videos, Ji et al. [16] showed that 3D CNNs sometimes improved performance over other methods including SVMs and 2D CNNs. Together, these studies suggest that better crop type classification performance can be achieved with deep learning methods that combine both recurrent and convolutional layers to model patterns in the spectral, spatial, and temporal dimensions of remote sensing observations. Many prior crop type classification studies have not assessed model generalization to future years that may exhibit growing season timelines that differ from the years used in training. Additionally, most prior studies rely on observations spanning the complete growing season, and thus do not provide within-season predictions of crop types. These limitations are barriers to operational adoption of crop type classification techniques. In this study, we present an approach for within-season crop type classification that is robust to shifts in growing season timelines by dynamically selecting input observations based on crop phenology. To assess the performance of this approach in future years that exhibit domain shift with respect to the training data, we evaluated our method for the 2019 growing season in the midwestern US in which planting was delayed by as much as 1-2 months.
We focused our study on a region in northern Illinois that covers approximately 12,056 km 2 including Knox and Stark counties as well as most of Peoria, Henry, and Bureau counties. We chose this region because northern Illinois was one of the regions most affected by extreme weather and floods in early 2019, and Illinois is the second-largest producer of corn in the United States [23] . Corn and soybeans are the dominant crops grown in this region, while forests, grassland, and urban development cover most area not covered by corn and soybean fields.
We used satellite observations from the Harmonized Landsat-8 and Sentinel-2 (HLS) and Sentinel-1 datasets. HLS adjusts Sentinel-2 surface reflectance products to have the same bandpasses and spatial resolution (30 m/pixel) as Landsat-8. This allows the Landsat-8 and Sentinel-2 satellites to be treated as a virtual constellation with 30 m/pixel spatial resolution and 2-3 days revisit time [8] . Sentinel-1 acquires C-band synthetic aperture radar (SAR) observations every 6 days at the equator. We used observations acquired between 2017-2019 since the second Sentinel-2 satellite was not launched until 2017. For Sentinel-1, we used VV and VH polarisation bands (speckle-filtered and normalized by incidence angle). For HLS, we used the blue, green, red, near infrared (NIR), and short-wave infrared (SWIR) bands (Landsat-8 bands B01-B06 and Sentinel-2 bands B02, 03, 04, 08A, 11, and 12) as well as two spectral indices: normalized difference water index (NDWI) and land surface water index (LSWI). We computed NDWI and LSWI using 41 Figure 1 : Inset: study area in northern Illinois covering 12,056 km 2 including Knox, Stark, Peoria, Henry, and Bureau counties (tile 16TBL in HLS reference system).
the following equations:
We used NDWI and LSWI instead of more commonly used indicese.g., normalized difference vegetation index (NDVI) or green chlorophyll vegetation index (GCVI)-because water content inside the leaf has been shown to an important factor for discriminating between corn and soybeans [5, 9] , which is captured by the SWIR bands included in the calculations for NDWI and LSWI. In our experiments we found that classification performance decreased when NDVI and GCVI were included in the input. We used the quality assessment (QA) layer provided in the HLS dataset to remove pixels with clouds and cloud shadows (bit numbers 0-3). The HLS tiles also contain no-data pixels (NaNs) in regions not covered by the satellite track for a particular date. To fill in no-data pixels and mitigate outliers in the time series data, we smoothed the time series for each pixel and each band using the Savitzky-Golay method [7] . We used the Cropland Data Layer (CDL) produced by the United States Department of Agriculture (USDA) National Agricultural Statistics Service (NASS) [2] as a surrogate for ground truth data, as has been done in prior studies (e.g., [6, 43] ). According to the CDL, approximately 60% of our study area is covered by corn or soybean fields, and most fields follow a corn-soybean rotation. Labels sourced from the CDL cannot be considered ground truth since they are predicted using a combination of ground truth labels and a decision tree model. However, the CDL typically has high accuracy for corn and soybeans in midwestern states [3] , with user's and producer's accuracies between 90-98% for corn and soybeans in Illinois estimated for CDL 2017-2019 [24] . Lark et al. [20] recommends post-processing steps to reduce the number of errors in the CDL labels, including aggregating classes and filtering out withinfield speckle and other heterogeneity. We aggregated all classes except corn and soybeans into an "other" class. We applied a 3 × 3 homogeneity filter to identify non-homogeneous regions (caused by within-field speckle) to be masked out during training.
Time series observations of crops from satellite images may have different temporal patterns from year to year depending on the planting timeline followed each year, which may change due to extreme weather or other climate factors. While methods that use the time series from the full growing season have been shown to have good performance for discriminating crop types, these methods will likely have poor performance for seasons in which planting timelines are different than in the training seasons (i.e., domain shift). Furthermore, most methods do not enable predictions of crop types within the growing season, when uncertainty around crop production is most uncertain. We propose a new approach for crop type classification that addresses domain shift by normalizing inputs by crop growth stage. In each pixel, we automatically select observations at key growth stages rather than using observations from the same dates for all pixels. These key observations are then fed to a deep neural network that combines recurrent and convolutional layers to capture spectral, spatial, and temporal relationships in the input observations.
While specific growth stages vary between crop types, all crops generally contain three key growth stages: greenup, during which crops emerge from the ground and undergo early vegetation; peak growth, during which crops have reached their maximum growth and leaf area index (LAI), and begin to ripen; and senescence, during which crops are drying down before they are harvested. For many crops including corn and soybeans, these growth stages are apparent in the time series of vegetation indices like NDVI (e.g., [1, 44] ), though regional planting timelines may vary from year to year due to weather or other climate-related factors that influence soil conditions or temperature [12] . For this reason, it is important that observations are compared at the same growth stage, which may not be at the same time [34] . Observations of the same crop on the same day but in different years will not be comparable if the planting timelines each year were different, resulting in domain shift. Figure 2 shows the distribution of NDVI values in corn, soybean, and other pixels as well as true color images on the same day of year (DOY) in 2017, 2018, and 2019. In 2019, the distribution of NDVI values at the same DOY is significantly lower for both corn and soybeans due to planting delays caused by extreme flooding, while the other class distribution showed only minor differences (likely due to minor crops that were also delayed).
Since observations at the same growth stage should be comparable regardless of planting timeline, our approach is to distinguish crop types using satellite observations acquired at key growth stages rather than using fixed dates or the full growing season. For each pixel in our dataset, we computed the approximate date of greenup, peak growth, and senescence using the NDVI time series. We defined the greenup DOY to be the day after which the slope of the NDVI time series curve is highest (i.e., increasing fastest) within the first half of the time series. We defined the DOY of senescence to be the day (following the greenup DOY) after which the slope of the NDVI time series curve is the lowest (i.e., the most negative, thus decreasing fastest). We defined the DOY of peak growth to be the day at which NDVI reaches its maximum value between the greenup and senescence DOYs. Figure 3 shows an example NDVI time series with greenup, peak growth, and senescence DOYs.
Prior studies detailed in Section 2 suggest that neural network architectures that use both convolutional and recurrent layers to capture spatial and temporal patterns in data (respectively), in addition to using observations from multiple spectral bands, should have the best performance for crop type classification compared to alternate methods. Our model architecture (shown in Figure 3 ) consists of two branches: a 2D CNN branch to capture spatial and/or textural patterns, and an LSTM branch to capture temporal patterns. The CNN and LSTM branches are joined by a dense layer before the final softmax prediction layer. Since the number of growing days is a distinguishing factor between corn and soybeans, we included a third input to the dense layer joining the CNN and LSTM branches that is the scalar number of days elapsed between the senescence DOY and greenup DOY, i.e.: δ = DOY senescence − DOY greenup . Recurrent branch. The first layer in the recurrent branch is a 1D convolution layer (8 filters, 1 × 3 kernels) followed by four LSTM layers with 16 units in all layers. All LSTM layers are densely connected, meaning the output of all previous LSTM layers is included in the input in addition to the output of the prior layer. Each LSTM layer is followed by a batch normalization [14] layer and dropout [39] layer (randomly dropping 30% of units). We used the rectified linear unit (ReLU) activation function for all layers. The input to the LSTM branch consists of 6 HLS spectral reflectance values, 2 HLS spectral index values, and 2 SAR polarisation values (see Section 3 and Figure 3 ) for each pixel acquired on the greenup, peak, and senescence DOYs-i.e., x LSTM i ∈ R M ×D where M = 10 is the number of channels and D = 3 is the number of timesteps.
Convolutional branch. The convolutional branch consists of two 2D convolution layers, each having 64 filters and 3 × 3 kernels with 1-pixel strides. We found that fewer or more than two convolutional layers did not improve performance for the validation dataset. We used ReLU activation for both layers, and each layer is followed by a batch normalization layer (dropout did not improve model generalization in this branch). We included a 2 × 2 max pooling layer after the two convolution layers. The input to the CNN branch consists of 5 × 5 image patches in which the central pixel is the pixel to be classified for all 10 channels and 3 timesteps. We also experimented with larger patch sizes (odd dimensions between 7 × 7-13 × 13), but found this did not improve validation accuracy. The multispectral patches for each timestep are concatenated to produce a 3D input tensor: x CNN We implemented our CNN-LSTM architecture using the Keras deep learning library with the TensorFlow backend. The code for our model and experiments (detailed in Section 5) is available at https://github.com/nasaharvest/croptype-mapping-gsn.
Our goal for all experiments in this study is to assess model generalization to observations acquired in future years that were not part of the training dataset, in order to assess the resilience of each method to time domain shifts. We used observations from 2017 and 2018 in the training and validation datasets and from 2019 in the test dataset, since 2017 and 2018 follow typical planting timelines for corn and soybeans in the US midwest while 2019 saw delayed planting (discussed in Section 1). We partitioned our study region (shown in Figure 1 ) into quadrants and used pixels in the northwest, northeast, and southeast quadrants for training and the southwest quadrant for validation. In all experiments, we used the validation dataset to assess when learning plateaued and ensure the model did not overfit to the training dataset. We conducted three experiments to evaluate the performance for crop type classification in the beginning, middle, and end (before harvest) of the growing season.
In the late-season classification experiment, classifier inputs contain observations acquired at all three growth stages: greenup, peak growth, and senescence. The mean senescence DOY for 2017-2019 is 278 ± 28. This is around the time that crops are harvested in this region (September-November) and production estimates remain uncertain, thus our model predictions would enable estimates several months before the CDL is made available publicly by the USDA (typically February the following year). To assess the contribution of each component in our model, we compared the performance of our model (denoted CNN-LSTM δ ) with the CNN branch only, the LSTM branch only, and the CNN-LSTM without the scalar δ input. To assess the benefit of normalizing inputs by growth stage rather than using fixed DOYs, we evaluated model performance using inputs from fixed DOYs (except CNN-LSTM δ since δ would be constant across all inputs). For the first date we used DOY 153, 149, and 152 for 2017, 2018, and 2019 respectively since greenup typically occurs around late May/early June; for the second date we used DOY 210, 209, and 211 since peak growth is typically reached around late July; and for the third date we used 274, 271, and 274 since senescence typically occurs around late September. Since random forests are widely used for crop type classification in prior work [26] , we also compared our method to a random forest with 900 estimators (chosen based on highest accuracy from 5-fold cross validation for estimators between 100-1000 with increments of 100). The input to the random forest included the 6 optical bands, 2 SAR bands, 2 spectral index values, and δ . We implemented the random forest using Scikit-learn [29] in python.
We trained all models for this and subsequent experiments for 25 epochs (after which validation accuracy plateaued) using minibatch sizes of 4096, cross entropy loss, and the Adam optimizer Table 1 : Performance metrics (overall accuracy (OA), per-class precision/user's accuracy (UA) and recall/producer's accuracy (PA)) for each crop type classification method using input observations from fixed and growth stage-normalized (GS-Norm) dates. Bold indicates best performance across all methods. with default settings [18] . To select training samples, we randomly sampled 520,000 pixels from the available training pixels (i.e., those not masked out by the homogeneity filter) and 10,000 from the available validation pixels. To mitigate class imbalance, we subsampled majority classes to have the same number of samples as the minority class (96,459 training and 5,517 validation samples per class). We used the same random sample and seed for all experiments. We observed empirically that validation accuracy decreased when a substantially larger or smaller number of training examples was used. We report the performance metrics (overall accuracy as well as per-class user's accuracy/precision and producer's accuracy/recall) for the five models using both growth-stage normalized inputs and fixed DOYs in Table 1 . We used the Sieve tool in QGIS to reduce within-field speckle in the CDL labels before computing performance metrics. In Figure 4 , we illustrate the qualitative performance in the predicted classification maps for our study area compared to the CDL reference labels. Table 1 shows that the growth-stage normalized CNN-LSTM δ model has the highest overall accuracy for the 2019 test data. Table  1 also shows that for all models, growth-stage normalized inputs give better performance than fixed DOYs. The CNN-LSTM architecture consists of two branches-an LSTM and 2D-CNN branch-with outputs joined by a dense (fully-connected) layer. The only difference between this and CNN-LSTM δ is the addition of a third input to this dense layer, which is the scalar difference between the senescence and greenup DOYs, intended to capture the number of growing days. Table 1 shows that including δ increases accuracy by 2.6%. The performance of the CNN branch-only model performs much better (by 5.8%) than the LSTM branch-only model, which Table 2 : Performance metrics (overall accuracy, per-class precision/user's accuracy (UA) and recall/producer's accuracy (PA)) for each crop type classification method for early-season and mid-season classification experiments. Bold indicates best performance across mid-season classification and early-season classification methods. suggests that spatial information is more informative for discriminating crop types than temporal information in this experiment. While prior work found recurrent layers to have greater influence on performance than convolutional layers (e.g., [4, 10] ), the difference is likely due to the shorter time series (3 timesteps) used in this study. The random forest had lower overall accuracy than the CNN-LSTM δ , but higher overall accuracy than the other neural network models. This is likely because the random forest includes the scalar δ input while the CNN-LSTM, CNN, and LSTM models do not, but does demonstrate that a random forest can achieve comparable performance to deep learning methods for this problem. Table  1 also shows that the difference between the random forest with growth-stage normalized and fixed DOY inputs is much lower than for the neural network models. This suggests that random forests may extrapolate better to changes in planting timelines than neural network approaches, since the random forest does not explicitly model temporal or spatial patterns that rely on the ordering of inputs.
Since all models in the late-season classification experiments described in the previous section performed better with growth stagenormalized inputs, we did not report results using fixed DOYs in the two remaining experiments for early-season and mid-season classification. We used the same experimental setup as for the lateseason classification, but only included greenup and peak growth stage inputs in the mid-season classification and only greenup for the early-season classification experiments. For the CNN-LSTM δ model, we computed δ as the difference between greenup and peak DOYs instead of greenup and senescence as in the late-season model. We only evaluated the CNN model and the random forest for earlyseason classification since there is no temporal information to be captured by recurrent layers. Table 2 gives the performance metrics for each model in the mid-season and early-season classification experiments. While the random forest results in Table 1 were achieved including δ in the model input, we found that the random forest had significantly higher accuracy for the mid-season classification experiment when δ was not included in the input (32.8% with δ , 76.5% without). Thus the random forest results in Table 2 were achieved without the δ input. As expected, the overall accuracy for all models decreased when the senescence DOY was excluded from the input in the mid-season classification experiment. This drop in accuracy is primarily due to low accuracy for the corn and soybean classes, as we observed that the recall for the other class remained high regardless of the number of timesteps. The CNN had the highest overall accuracy (82.8%) for mid-season classification, emphasizing that recurrent layers are less effective given fewer temporal observations. The mean peak growth DOY in our experiments is 216 ± 27, so this classification performance could be achieved between July-September within the growing season.
For the early-season classification experiment, the overall accuracy is significantly lower when only the greenup timestep is included in the input to the CNN (27.3% lower than in mid-season). This underscores the conclusion from prior studies and this study that crop type classification performance improves significantly when observations from multiple timesteps are used (even if only a few dates). The random forest accuracy is also lower in early-season than mid-season (7.5% lower), but this difference is much lower than for the CNN. We also found that with only one timestep, the random forest performs better than the CNN. Thus, while the neural network approaches are more effective when additional growth stages are included later in the season, the random forest has the best performance given limited data in the beginning of the season. This could also be a result of the complexity of the neural network models, which is more suitable for the higher-dimension inputs in the later-season experiments. The mean greenup DOY in our experiments is 144 ± 37, so early-season predictions could be made using the random forest with 69.0% accuracy between April-June near the start of the growing season.
Most prior work on crop type classification has focused on postseason classification and have not addressed the effect of domain shift on model performance in future years not included in the training dataset. To be used operationally, models must be resilient to domain shifts that may occur in future seasons-e.g., due to delayed planting-as the global climate changes and extreme weather events become more frequent. Operational methods also can't rely on collecting more data each year to capture these shifts in training examples, e.g., due to the COVID pandemic or in countries where security is compromised. Rather than using observations from a fixed time period or set of dates, our approach normalizes inputs to be captured at the same key crop growth to ensure comparability between observations acquired in different years regardless of the planting timeline. Our experiments in Section 5 showed the best within-season classification accuracy of 85.4% was achieved using our proposed CNN-LSTM δ model with growth-stage normalized inputs late in the season prior to harvest (September-November), and that overall accuracy of 82.8% could be achieved in mid-season (July-September) and 69% at the beginning of the growing season (April-June) using the growth-stage normalized CNN and random forest respectively. In future work we will determine if this remains true when the models are applied to larger, more geographically diverse areas in the US Corn Belt states.
Our experiments confirmed prior findings that neural network architectures that used both recurrent and convolutional layers with multispectral, multi-temporal observations gave the best crop type classification performance (e.g., [4, 10] ). Our late-season classification model included inputs from only three timesteps; in future work we will explore using more timesteps between the greenup and senescence DOYs to make better use of the recurrent layers. Although one of the main benefits of neural networks is that useful features can be automatically learned from data without the need for domain-specific feature engineering, we found that introducing some domain-specific information in the form of growth stage-normalized inputs enabled all models (including random forests) to be more robust to patterns in future (test) seasons not captured in the training data. Consistent with prior studies (e.g., [11, 17, 19, 26, 32, 37] ), we also found that including SAR inputs from Sentinel-1 in addition to optical inputs significantly improved classification accuracy; in future work, we will explore additional SAR parameters in addition to VV and VH polarization that can be extracted from SAR data. Finally, we plan to investigate aggregating pixel-level observations or predictions at the field-level, which has been shown to improve classification performance in prior work (e.g., [5, 6] ).
Crop type classification is an important technique for providing information about planted area and enabling estimates of crop condition and forecasted yield, especially within the growing season when uncertainties around these quantities are highest. As the global climate changes and extreme weather events become more frequent, crop type classification methods (for both within-season and post-season classification) must be resilient to domain shifts (e.g., due to changes in planting timelines) that can adversely affect predictive performance. Models that provide reliable predictions without the need to collect new training data each season is especially important when ground surveys are difficult or impossible, e.g., due to conflict-related insecurity or travel restrictions imposed from the COVID-19 pandemic.
We proposed an approach for within-season crop type classification that uses remote sensing inputs normalized by crop growth stage to ensure compatibility across seasons even when domain shifts may occur. We used a neural network leveraging both convolutional and recurrent layers to predict whether pixels belong to the corn, soybean, or other class. We evaluated this method for the 2019 growing season in the midwestern United States, during which planting was delayed by as much as 1-2 months due to extreme weather that caused record flooding. We showed the highest classification accuracy of 85.4% could be achieved near harvest between September-November, 82.8% accuracy mid-season between July-September, and 69.9% between May-July in the beginning of the growing season. Our study was limited to a region spanning several counties in northern Illinois, and our next step is to scale our approach to all of the US Corn Belt states. We will assess operational use of this approach during the 2020 growing season to provide information to help mitigate uncertainty in agricultural markets caused by the COVID-19 pandemic.