Combining GEDI and Sentinel-2 for wall-to-wall mapping of tall and short crops

To evaluate the potential of GEDI to distinguish between tall and short crops, we considered three regions of the world: Jilin in China, Grand Est in France, and Iowa in the United States (figure 2). These regions are major agricultural production areas located on three separate continents. They each contain a mix of tall and short crops and have accurate, up-to-date field-scale crop type maps that are publicly available.

Zoom In Zoom Out Reset image size

Jilin Province is located in Northeast China and is one of the most important maize-producing provinces in China. It spans the mid-latitudes from 40.8^∘N–46.3^∘N and 121.7^∘E–131.3^∘E and has a humid continental climate. Other major crops grown in the area are soybeans and rice. Because early frost usually appears in September and early October, fast maturing maize varieties are cultivated. Maize is typically planted in April and harvested in September, with an average maize cycle duration of 150 days.

The Grand Est administrative region in northeastern France cultivates a wide variety of crops including wheat, barley, maize, alfalfa, sugar beets, legumes, and oilseeds. Twenty-two percent of French sugar production, 21% of rapeseed production, 13% of wheat production, and 13% of maize production come from this region (Le Service statistique ministériel de l'agriculture 2019). Maize is generally planted between April and May, and harvested between September and November. The region spans 47.4^∘N–50.2^∘N and 3.4^∘E–8.2^∘E and has a climate that varies from oceanic in the west to humid continental in the east. While the administrative region Nouvelle-Aquitaine is the largest producer of maize in France (31%), Nouvelle-Aquitaine also produces significant quantities of sunflower, which is also a tall plant. To focus on evaluating GEDI's ability to distinguish maize, we instead conducted experiments in Grand Est, which is France's second-largest producer of maize. We elaborate on the application of GEDI in regions with more than one tall crop in the Discussion.

Iowa, a state located in the Midwestern region of the United States, is in the heart of the U.S. Corn Belt and is the country's largest producer of maize. Maize and soybean are the two primary crops cultivated, comprising well over 95% of total cropped area (figure A1). Maize in Iowa is planted from late April to May and harvested from late September to early November. Located between 40.4^∘N–43.5^∘N and 90.1^∘W–96.6^∘W, Iowa also experiences a humid continental climate with cold winters and hot summers.

While the three study regions are thousands of miles apart, they share some similarities: they are all in the Northern Hemisphere, have humid continental climates, grow predominantly rainfed agriculture, and sow maize around April and harvest it between September and November. Management practices, however, do differ; of the three study regions, Jilin has the smallest fields on average, with the median field size between 0.64 ha and 2.56 ha according to the GeoWiki database (Lesiv et al 2019). Fields in Grand Est are 4.5 ha on average (Agence de Services et de Paiement 2019b), while fields in Iowa are largest with an average area of 33 ha (Yan and Roy 2016). Differences in agricultural practices across regions may affect model transfer if the practices result in the same crop type having different GEDI waveforms. Smaller field sizes also cause more GEDI shots to land on field borders, which are likely to contain more than one crop type or a mix of crop and non-crop land covers.

The GEDI instrument is the first spaceborne lidar instrument specifically optimized to measure vegetation structure. By firing a laser at 25 meter spots (termed 'footprints') on the Earth's surface and observing the return of the laser pulse, the instrument is able to measure the vertical distribution of vegetation at each spot. GEDI collects data globally (between 51.6^∘N and 51.6^∘S latitudes) at the highest resolution and densest sampling of any lidar instrument in orbit to date (Dubayah et al 2020). The raw GEDI waveforms collected undergo several processing steps to retrieve a variety of metrics, and the derived products are saved as both footprint and gridded datasets.

For our analysis, we used the Level 2A Elevation and Height Metrics Data (L2A), which includes footprint-level elevation and relative height (RH) metrics. RH metrics represent the height (in meters) at which a percentile of the laser's energy is returned relative to the ground. For example, $\textrm{RH}50 = 20$ m means that 50% of the laser's energy was returned by objects up to 20 m above the ground. The ground position is determined based on the center of the lowest mode of the returned waveform (Hofton et al 2000). RH metrics are saved at 1% intervals, so each shot contains 101 values representing RH at 0%–100%. These footprint data are geolocated with a mean positional error of 10.3 m.

We downloaded the GEDI L2A version 2 (GEDI02_A v002) data from July to September 2019 that intersects our study regions through NASA's Earthdata Search website. GEDI L2A data come with a series of flags and properties to help the user filter for data of quality appropriate for the specific application. For this study, we omitted shots with a quality flag value of zero, which indicates poor quality, and a non-zero degrade flag, which indicates poor geolocation. We note that, unlike RH values observed for forests, RH values used here were commonly below zero. This happens because waveforms from agricultural areas often have only one mode, and the GEDI algorithms define RH relative to the center of the lowest mode. We also filtered out shots with RH100 greater than 10 m, as the field crops in the study areas do not grow that tall. We also dropped a full orbit for September 25 in Jilin, China because of abnormally high RH100 values. The outliers removal only dropped a small percentage of points (1.5% in Iowa, 3.6% in Jilin, and 3.5% in Grand Est). A map of the shots left in each region after cleaning the dataset and filtering for cropland are shown in figure 2, and the counts of shot numbers for each crop type are summarized in figure A1.

Consecutive RH metrics are highly correlated with each other. We therefore sampled a metric every 10% to reduce the number of features used from 101 to 11. The accuracy of random forest classifiers trained on a subset of 11 RH metrics is only slightly lower than one trained on all 101 RH metrics for all three study regions. This suggests that the information lost from reducing feature dimensionality had little impact on the ability to distinguish crop types.

Current state-of-the-art crop type maps use time series from passive optical remote sensing as input features for classification (Defourny et al 2019, USDA-NASS 2020). To compare GEDI features to optical features for crop type classification, we extracted Sentinel-2 (S2) time series at each GEDI shot location. We also extracted S2 time series for the entirety of the study areas in order to demonstrate how GEDI can be used to create labels for wall-to-wall crop type mapping. All optical imagery was processed using the Google Earth Engine (GEE) platform.

The Sentinel-2A/B (S2) satellites acquire images with a spatial resolution of 10 m (Blue, Green, Red, and NIR bands) and 20 m (Red Edge 1, Red Edge 2, Red Edge 3, Red Edge 4, SWIR1, and SWIR2 bands), and together they provide images at a 5-day interval. The spatial resolution of 10 to 20 m is sufficient to resolve individual fields in the three study areas.

We used S2 surface reflectance data (Level-2A) present in GEE and filtered out clouds using the S2 Cloud Probability dataset provided by SentinelHub in GEE. To capture crop phenology, we used Sentinel-2 imagery from 1 January to 31 December 2019, using the same time window across the three regions. In our study areas, this time window encompasses a single growing season for the majority of crop types.

Features were extracted from S2 time series by fitting harmonic regressions to all cloud-free observations in 2019. The harmonic regression, also known as the discrete Fourier transform, decomposes signals over time into their constituent frequencies (cosines and sines). It has been shown to perform well at land cover and crop type classification (Moody and Johnson 2001, Jakubauskas et al 2002). Other options for feature extraction include monthly or seasonal medians and the double logistic regression. In our previous work in the U.S. Midwest (Wang et al 2019), we found that harmonic coefficients classified crop types better than seasonal median features. Furthermore, unlike the double logistic regression, the harmonic regression is easy to deploy over entire states and provinces in GEE using the built-in linear regression function.

For each spectral band or vegetation index f(t), the harmonic regression takes the form

$$\begin{equation} f(t) = c + \sum_{k = 1}^{n} \left[ a_{k} \cos (2~\pi \omega k t) + b_{k} \sin(2~\pi \omega k t) \right] \end{equation} \tag{ 1 }$$

where a_k are cosine coefficients, b_k are sine coefficients, and c is the intercept term. The independent variable t represents the time an image is taken within a year expressed as a fraction between 0 (1 January) and 1 (31 December). The number of harmonic terms n and the periodicity of the harmonic basis controlled by ω are hyperparameters of the regression. We used a second order harmonic (n = 2) with ω = 1.5, shown in previous work (Wang et al 2019) to result in good features for crop type classification in the U.S. Midwest. Since seasonality is similar in Jilin and Grand Est, the same values of ω and n should also fit time series in the other two regions well. Model transfer also requires that the features be derived from the same regression. The harmonic regression above yields a total of five features per band or vegetation index, resulting in 20 harmonic coefficients total.

We computed harmonic coefficients for three bands and one vegetation index: NIR, SWIR1, SWIR2, and GCVI. GCVI is the green chlorophyll vegetation index (Gitelson et al 2005) computed as

$$\begin{equation*} \textrm{GCVI} = \textrm{NIR} / \textrm{Green} - 1. \end{equation*}$$

Unlike the commonly-used NDVI, GCVI does not saturate at high values of leaf area. These four bands and VI were chosen in prior work (Wang et al 2019), which found that crop type classification using NIR, SWIR1, SWIR2, and GCVI performed nearly as well as classification using all optical bands and a variety of other VIs.

We used high-accuracy crop type maps in Jilin, China, Grand Est, France and Iowa, USA to filter out non-crop areas, train maize classifiers, and evaluate each classifier's performance. In each region, the corresponding map's value at each GEDI shot footprint centroid location or S2 pixel location was used as the ground truth for crop type. Note that GEDI footprints have a 12.5 m radius, so it is possible for a GEDI shot to span multiple crop type map pixels and have mixed crop type labels (figure 3). The data products available in each region for the year 2019 are described below.

Zoom In Zoom Out Reset image size

2.4.1. You et al (2021) Northeast China crop type map

2.4.2. Registre parcellaire graphique (RPG)

2.4.3. USDA CDL

We used a random forest classifier to classify maize vs. non-maize in all experiments. Random forests (Ho 1995, Breiman 2001) are an ensemble machine learning method comprised of many decision trees in aggregate. Each decision tree is trained on a bootstrapped version of the training set and a random subset of features to reduce the correlation of predictions across decision trees and improve performance when those predictions are averaged. Random forests are commonly used in crop type classification (Defourny et al 2019, Jin et al 2019) and other Earth observation tasks due to their high accuracy and computational efficiency.

We used the RandomForestClassfier implemented in Python's scikit-learn package. We kept the default parameters, with the exception of raising n_estimators from 10 to 100 to reduce prediction variance.

We discretized each study region into $0.5^{\circ} \times 0.5^{\circ}$ grid cells. All GEDI shots in each grid cell were placed entirely in either the training set or test set, with 80% of grid cells in training and 20% in test (figure A2). Splitting GEDI shots along grid cells reduces spatial correlation across the training and test sets, thereby preventing classification metrics from being artificially inflated. For each result reported in this paper, we ran the classifier 11 times using different training and test splits, and we report the mean and standard deviation of accuracy over all runs. We chose an odd number of runs (11 instead of, say, 10) to make it easier to select the run with median accuracy and visualize its confusion matrix.

Our first experiments test how well GEDI waveforms can distinguish maize from non-maize crops based on maize being a significantly taller crop. As GEDI was designed to monitor forests, it is unknown whether the instrument would be able to resolve height differences between crop types at all.

We trained a random forest classifier for each study region using GEDI RH metrics ('GEDI Local') as features, and compared this model to a baseline random forest classifier trained on S2 harmonic coefficients ('S2 Local'). Both models were tasked with distinguishing maize samples from non-maize samples, where the ground truth for crop type was provided by the datasets described in section 2.4. The S2 Local model is representative of current state-of-the-art crop type classifiers that use optical imagery to predict crop types. It provides a reference for the GEDI Local model as well as models that are transferred across regions.

The number of samples and locations used for training and testing the GEDI and S2 models were identical. Although S2 provides wall-to-wall imagery while GEDI only observes a small subset of Earth's surface, we limited S2 samples to GEDI shot locations. By controlling for sample size and location, we directly compare the two sets of features.

The timing of GEDI observations is important, as maize will be most distinguishable from other crops when their height difference is the greatest. To test the sensitivity of classification performance to growing season timing, we compared the performance of GEDI Local models trained on July-only shots, August-only shots, September-only shots, and shots from all three months.

After classifying maize vs. non-maize crops within each region, we tested model transfer across regions. For each GEDI Local and S2 Local classifier trained in the U.S., China, or France, we applied the classifier to the test sets of the other two regions to separate maize from non-maize. We refer to these models as 'GEDI Transfer' and 'S2 Transfer'. The models were not shown any additional data from the new regions. High classification accuracy in a new region would indicate that the model's features generalize across space and few if any labels are needed from the new region to classify maize. Conversely, low classification accuracy would mean that the learned relationship between features and crop types holds true only locally, and labeled data is needed in the new region to learn new classification boundaries.

We compared the GEDI Transfer models to their Local counterparts to see how well GEDI RH metrics generalize across geography. To understand whether growing season timing affects model transfer, we repeated this analysis for each GEDI model trained on July-only shots, August-only shots, September-only shots, and shots from all three months.

Even if GEDI features transfer perfectly across geography—i.e. maize is always identifiable in GEDI waveforms no matter where on Earth one looks—GEDI only samples 4% of the land surface and cannot alone generate a wall-to-wall crop type map. Achieving a continuous map in space requires GEDI-based approaches to be combined with wall-to-wall imagery like that provided by S2.

To create wall-to-wall maize maps using GEDI and S2 imagery, we used methods similar to those employed previously to calibrate local maps of forest height with GEDI and Landsat (Healey et al 2020). For each pair of regions, we first followed the steps in sections 3.2 and 3.3 to obtain GEDI Transfer models. For example, in China we had one GEDI Transfer model trained in the U.S. and another one trained in France. Then, for each GEDI Transfer model, we treated its predictions at GEDI shot locations as labels in the new region and trained another model to reproduce those predictions using local S2 harmonics as features. We call this second model 'GEDI-S2 Transfer'. To continue the example, we treated the predictions of the U.S.-trained (or France-trained) GEDI Transfer model at GEDI shots in China as true crop type labels, and then trained a classifier on S2 features in China to predict those 'labels'.

By applying the GEDI-S2 Transfer model to all cropland S2 pixels in a new region, we produced a wall-to-wall 10 m spatial resolution maize map without the need for local labels. Figure 4 presents a graphical explanation of the GEDI-S2 Transfer approach with transfer from the U.S. to China as a concrete example.

Zoom In Zoom Out Reset image size

We compared GEDI-S2 Transfer to two other models. The first is the S2 Local model, which provides an upper bound for how well S2 harmonic features can classify maize when trained on in-region ground truth. The gap between GEDI-S2 Transfer and S2 Local reflects the accuracy of GEDI Transfer's predictions relative to ground truth. The second benchmark is the S2 Transfer model, which shows how well a model trained on S2 harmonics in one region fares when applied to other regions. Differences between S2 Transfer and GEDI-S2 Transfer reveal how robust GEDI features are across space compared to optical features.

The median harmonics of GCVI from S2 and median RH energy curves from GEDI are shown for the top three crops in each region in figure 5. In each region, the S2 maize profile reaches peak greenness in August and is generally distinguishable from the S2 profiles of other crops because of different timing and magnitude of the greenness peaks. For the GEDI energy profiles, roughly 50% of the energy returned comes from negative RH values, which as mentioned previously arises from the GEDI algorithm's definition of ground elevation (section 2.2). Thus, we emphasize that the values of RH should not be interpreted as physically meaningful; for instance RH100 does not correspond to the physical crop height. Nonetheless, the GEDI curves exhibit a clear separation between maize—the tallest crop—and the other crops. This difference is especially apparent at the extremes of RH curves, shown as insets in figure 5. Specifically, maize has lower values for RH0 to RH50 and higher values for RH70 to RH100 compared to other crops. These differences are also apparent in the example profiles shown in figure 3.

Zoom In Zoom Out Reset image size

Whereas figure 5 compares different crops within a region, figure 6 displays the median features for maize from different regions on the same plot. This comparison is especially relevant for the question of how well a model is likely to transfer across regions. Ideally, features would be similar across regions in order to use a model trained in one region on another. For the S2 harmonics, clear differences emerge between the regions. Maize in the U.S. generally has a steeper increase in GCVI during June and July and a higher peak in August compared to maize in the other regions. Maize curves in France are slightly earlier than maize curves in China on average and have a substantially higher variance. The harmonic curves for the non-maize crops also differ considerably, both in timing and magnitude of the peak (figure 5). In contrast, the GEDI curves are remarkably consistent across regions (figure 6). This difference between S2 and GEDI indicates that the peak height of maize is a better preserved characteristic across regions than the timing of the maize growing season and the total crop biomass, both of which influence the S2 harmonics.

Zoom In Zoom Out Reset image size

The feature comparisons above suggest that GEDI features should be reliable for classifying tall crops (in this case, maize). To quantitatively evaluate its performance for local classification, we used these GEDI features to train a classifier for each region and compared it to a classifier trained on the S2 harmonics (section 3.2).

Since the task is to separate maize from other crops based on height, we considered various possible timings of the GEDI features. Test accuracies indicate that GEDI RH metrics can distinguish maize from other crops within each region with more than 79% accuracy in all three months tested (figure 7). The optimal timing of GEDI observations differs by region, with September for China, July for France, and August for the U.S. being the best times for classification, resulting in 88%, 85%, and 91% accuracy, respectively. August is generally a good month for GEDI observations in all regions, with performance in all regions above 83% for this month.

Zoom In Zoom Out Reset image size

As expected, the locally-trained S2 models do well in each region (93% accuracy in China, 95% in France, and 95% in the U.S.), and indeed perform better than the GEDI features. The gap between S2 local models (trained on the same locations as the July–September GEDI shots) and the best GEDI models is less than 5% for China and the U.S. and less than 10% in France.

To understand why GEDI model errors are larger, we show maps of GEDI misclassifications in appendix figure A3 and confusion matrices for representative median runs in appendix figure A4. From the maps, we see that a significant percent of errors occur at field borders, where the GEDI shot footprint contains multiple crop types or a mix of crop and non-crop classes. We also observe that the most difficult crop types for GEDI Local models to classify are soybeans in China and silage corn and sunflower in France. In China, this can be partly explained by the 84% precision for soybeans in the You et al (2021) map used for ground truth (section 2.4.1). In other words, up to half of the 'misclassification' of soybeans as maize in China could be correct. Since the 'ground truth' maps in China and the U.S. are created using optical satellites, the predictions of S2 Local models—both correct and incorrect—are likely to correlate with the ground truth more than those of GEDI Local models.

Local GEDI accuracies in France are overall lower than the other two regions. This appears to arise mainly from the fact that two kinds of maize are grown in France, maize for grain and for silage (see figure A1 for the crop distribution by region). Whereas grain maize is always grown to maturity and thus has a more reliable seasonality, silage maize is grown for biomass and can be sown and harvested at any point in the season. As a sensitivity test, we recalculate local GEDI accuracy in France after omitting the silage maize from the test set, finding that accuracies improve by 4% or more in all periods (see figure A5 in appendix). The improvement is largest in September, consistent with the notion that early harvest of silage maize is affecting the performance of GEDI features, since a harvested crop is no longer tall. As both the U.S. and China grow predominantly grain maize, the less predictable timing of silage maize is not an issue in those regions.

As noted in section 4.1, the consistency of GEDI features across regions suggests that models trained in one region can be reliably applied to new regions. A quantitative test of this proposition is shown in figure 8, which compares the performance of GEDI models trained using data from the local region (GEDI Local) to those trained in other regions (GEDI Transfer). Although the locally trained models are typically the best performers, the transferred models perform nearly as well and are occasionally indistinguishable from the locally trained models. For example, models trained in the U.S. or China both perform as well in France in August as one trained in France.

Zoom In Zoom Out Reset image size

Transferred GEDI models are only able to make predictions at GEDI shot locations. In order to extrapolate beyond these point locations, we used transferred GEDI predictions as labels to train a new model that takes local S2 harmonics as input (GEDI-S2 Transfer). GEDI-S2 Transfer test accuracies for the month of August are reported in figure 9 together with S2 Transfer and S2 Local accuracies for comparison. S2 Transfer shows how well state-of-the-art optical features perform out-of-region, while S2 Local provides an upper bound on how well optical features can separate maize when paired with ample local ground truth.

Zoom In Zoom Out Reset image size

The results in figure 9 show first that harmonic features, while good at distinguishing crop types locally, transfer poorly across study regions. The average accuracy of an S2 Transfer classifier is 64%. In the U.S., the S2 Local model achieves an accuracy of 94%, but S2 Transfer models trained in China and France only manage accuracies of 60% and 62%, respectively. This is consistent with previous work that found optical feature transfer-ability to deteriorate across geography (Wang et al 2019), as well as with the differences in crop phenology observed in figures 5 and 6. As the phenology and prevalence of crop types shift across regions, the optimal decision boundary for classifying maize versus non-maize with harmonic features also changes. S2 Transfer therefore results in many misclassified samples.

GEDI RH features, on the other hand, transfer much better across regions, and consequently the S2 model they supervise (GEDI-S2 Transfer) also performs much better than direct S2 Transfer. Figure 9 shows that GEDI-S2 Transfer accuracies exceed 82% for all cross-region pairs. For example, GEDI-S2 Transfer achieves 86% accuracy in the U.S. for models trained in China or France. While the S2 Local model has an 8% higher accuracy, the GEDI-S2 Transfer significantly outperforms S2 Transfer from China and France by 26% and 24%, respectively. GEDI-S2 Transfer and GEDI Transfer accuracies are about the same; in China they are the same, in the U.S. GEDI-S2 Transfer is 1% lower, and in France GEDI-S2 Transfer is 1% lower for the model trained in China and 1% higher for the model trained in the U.S.

Example crop type predictions for the three regions using the best GEDI-S2 Transfer model are shown in figure 10. The corresponding ground truth crop type maps are also shown for reference. For the most part, the GEDI-S2 Transfer predictions agree with the ground truth maps.

Zoom In Zoom Out Reset image size