Explainability of Sub-Field Level Crop Yield Prediction using Remote Sensing

We collected agricultural data from different regions in South America and Europe, as depicted in Figure 2. Specifically, we gathered soybean crop data from two countries: Argentina and Uruguay, involving 10 and 8 different farms, respectively. In Northern Germany, we used data from wheat and rapeseed crop fields, all originating from the same 6 farms. In Argentina, all soybean farms used in our study are situated in the northern regions, close to Uruguay, while in Germany, rapeseed and wheat crops are cultivated in the same farms, located in the northwest and close to the center of the country.

To gain a comprehensive understanding of the climate conditions and topography in these regions, we rely on the Global Agro-Ecological Zones (GAEZ) model documentation, specifically for insights into land, water ¹¹1Available from https://gaez-services.fao.org/apps/theme-1/, accessed 20 March 2024. and agro-climatic ²²2Available from https://gaez-services.fao.org/apps/theme-2/, accessed 20 March 2024. ressources [32]. All our study regions share a sub-tropic climate with moderately cool temperatures. Furthermore, the air is humid in North Argentina and Uruguay, sub-humid without significant soil or terrain constraints in German fields near the country’s center, and semi-arid with notable soil and terrain limitations in northern Germany. The total number of rain days (days with daily precipitation exceeding 1 mm) varies across these regions. In North Argentina and Uruguay, it ranges from 100 to 130 rain days, with annual precipitation averaging around 1900mm and 1000-1900mm, respectively. In Northern Germany, the number of rain days varies from 160 to 220, with annual precipitation ranging between 500 and 700mm. The length of the growing season is approximately 200 days in Germany and extends to 300 days in North Argentina and Uruguay. Note that the annual precipitation values and growing period lengths are averaged over the period 1981-2010. As for the topography of the study sites, the median altitude exhibits significant variations. In German fields, altitudes range from 40m in the northwest to around 250m near the country’s center. In Uruguay, altitudes span from 30m to 200m, while in Argentina, the altitude variation is relatively smaller, ranging from 30m to 120m.

In our study, we apply XAI to the ML-based crop yield prediction task on the subfield level. The target dataset consists of rasterized yield maps with a resolution of 10 meters, which we divided into four datasets, as outlined in Table I. It is important to note that a field can be counted more than once if its yield data was collected for multiple years. Uruguay provides the largest amount of soybean fields for our study, primarily due to the number of farms involved, followed by Argentina and Germany. The field area varies from one country to another, which explains the disproportionality between the numbers of fields and pixels. To clean the yield data, we implemented a yield thresholding at 15t/h on the yield maps, as recommended by experts in agriculture. We subsequently applied an outlier filter using the three-sigma rule, removing data points falling beyond the range defined by three standard deviations from the mean. Further cleaning included an inter-quartile-range filter, as detailed in [33].

To predict the yield, separate models are trained for each dataset. We supply the model mainly with satellite data, and conduct additional experiments adding more modalities to the network. More specifically, we use the entire spectral signal of Sentinel-2 Level-2A data and collect all scenes from seeding to harvesting. In terms of additional data sources (ADS), and following the same process as in [34], we collect information about soil condition, digital elevation maps (DEMs), and weather data. All maps were rasterized to match the yield resolution (10m). In Appendix A, a summary of the predictive features used is shown in Table VIII, along additional information about the processing of input and output data.

Feature attribution methods are commonly used interpretation techniques that, when combined with domain knowledge, provide explanations for individual predictions. An attribution method assigns sensitivity or contribution scores to the input features, with the absolute value of these scores quantifying their importance on the model’s prediction. In our paper, we will evaluate and employ such techniques to explain the yield prediction models.

Given the $i^{\text{th}}$ pixel sample $X^{i}\in\mathbb{R}^{B\times T}$, where $B$ and $T$ are the numbers of bands and time steps, respectively, let $a^{i}_{b,t}\in\mathbb{R}$ represent the attribution result for the feature $x^{i}_{b,t}$. The relative value of $a^{i}_{b,t}$ indicates the extent to which this feature contributed to increasing or decreasing the final predicted value $\hat{y}^{i}\in\mathbb{R}$.

As discussed in [41], the process of explaining a model through feature attribution methods requires a solid comprehension of the predictive features being employed. Conventional methodologies for crop modeling rely on the crafting of features by experts, with the aim of encapsulating the extensive array of variables recognized as direct yield drivers. However, when ML and DL models are used, the scenario changes. Raw information, acknowledged for its potential to serve as a proxy for crop growth and health, is directly supplied to the network with minimal to no alterations. In our work, we address this limitation of explaining ML models in two dimensions: inclusion of vegetation indices and aggregation of the temporal dimension.

With the aim of augmenting the interpretability of the spectral dimension, we employ well-established vegetation indices, frequently employed in the literature for crop classification and yield prediction. These vegetation indices merge two or more spectral bands and generate an index with agronomic relevance, rendering it more easily interpretable than the raw bands. The fundamental idea entails substituting the original S2 raw bands with pre-calculated indices during the training phase. Subsequently, we will evaluate the impact of this adjustment on the model’s performance. In our experiments, we will investigate the impact on accuracy when training the model using vegetation indices and compare it to training the model using raw inputs. If the results are similar, we will continue to use the model trained with vegetation indices for future experiments. However, if there is a significant decrease in performance, we will not further use the improvement in spectral interpretability. The full list of the indices used and their respective formulas are summarized in Appendix C.

To gain insights and derive explanations from crop yield models using time series, we can use different types of domain knowledge, with particular potential in the different phases that each crop goes through. To enhance the temporal interpretability of the input feature and derive explanations from it, we will exploit the principal growth stages of the BBCH scale for rapeseed and wheat [42], and a modified scale for soybean [43], as summarized in Appendix E. Concretely we first get the modified attributions $\bar{a}^{i}_{b,t}$ for each explained pixel sample $X^{i}$. We further define for each growth stage $t_{\text{gs}}$ a set $T_{\text{gs}}$ of the time steps it covers, then compute the total attribution as follows:

A major benefit of this approach is that it overcomes the misalignment of the seeding and harvesting months between the different fields in each dataset, which would have hindered the interpretability of the raw attributions on a higher level.

We first evaluate the performance of the yield prediction models, where one model is trained for each dataset and for each sampling and feature composition setting. The results, based on the evaluation metrics used, are summarized in Tables II and III for both the raw and monthly samplings, respectively.

We conduct the evaluation at both the subfield and field levels. For the former, we compute each metric for each pixel and then average over the entire dataset. In the latter case, we first average the predicted and target yield values over the pixels within each field. Then, we compute the metrics for each field and finally, we average them over the entire dataset.

In all our experiments, we consistently observe that the results at the field level outperform those at the subfield level. This disparity arises because subfield metrics measure the model’s ability to capture the in-field variability of yield values. In contrast, on the field level, we notice that averaging the individual predictions can more easily and accurately estimate the average yield per field.

When comparing the $\text{R}^{2}$ results across the different datasets, we note that the top scores at the field level are generally comparable across the datasets, except for GER-W, which consistently ranks last by a considerable margin. At the subfield level, the results achieved in ARG-S clearly outperform those in other datasets, regardless of the temporal samplings used, while GER-W is still not reaching comparable scores. These disparities may originate from the distinct agronomic characteristics of each crop. Therefore, it is possible that the input data in our experiments might be inherently insufficient to capture these unique characteristics and establish a strong connection with the target yield, resulting in less accurate predictions. Furthermore, satellite data availability can also vary from one country to another, influenced by factors such as geolocation and cloud coverage related to the climate in each region [44]. The quality of the yield data collected from different countries may also further impact the overall results.

However, when replacing the full spectrum of satellite data with vegetation indices, we observe different responses. In the monthly sampling, this change had a minor negative impact in URG-S and GER-R, but in contrast, it led to a substantial decrease in scores in ARG-S and GER-W datasets. Notably, in these two datasets, the $\text{R}^{2}$ score dropped by 10 and 34 percentage points (p.p), respectively, at the field level. When applying this same operation to raw time series data, a similar decrease in scores is observed in Germany’s datasets, with a 14 p.p drop in rapeseed fields and a 44 p.p drop in wheat fields at the field level. However, the use of vegetation indices yielded comparable results to the full satellite spectrum in Argentina and even slightly improved results at the field level in Uruguay. One reason for these results might be the need for complex spectral band interactions. Vegetation indices usually map only linear or low-complexity inter-band interactions, while neural networks are able to learn rich non-linear interactions.

Analyzing the impact of enriching the satellite data with additional modalities, as shown in Table III, reveals an overall improvement in the results across all datasets. This results also align with the literature [45]. The smallest increase is observed in the Argentina dataset, while the most significant impact is seen in the rapeseed dataset in Germany.

To identify a robust attribution method for the interpretation of the yield prediction models, we conducted a comparison of different techniques using the same dataset and model. Using the mean baseline vector, we compute both spectral and temporal importances, as illustrated in Figure 3. We observe a consistent pattern in the temporal attributions. Specifically, time steps before seeding and after harvesting months (approximately from time step 9 to 16) exhibited low importances, while the early growing months had a more significant influence, and the late stages were the most important. However, looking at the spectral importance, it is hard to decipher an agreement pattern across the feature attribution methods. This highlights the need for further evaluation tools to assess these methods better.

In Figure 4, we assess the stability of the methods and their robustness against both minor and significant input modifications. This is quantified using the infidelity and sensitivity scores. The boxplots show the distributions of both evaluation metrics over the attributions of all explained samples, when different attribution methods are used. Our goal is to identify a method for which the attributions yield low sensitivity and infidelity. Our observations reveal that most methods exhibit high sensitivity to minor input perturbations, with only four methods maintaining flat scores when subjected to significant noise in the input data, i.e. flat infidelity scores. Notably, SVS stands out as the sole method that remains robust against both perturbation types.

We illustrate this characteristic further in Figure 5 by examining the attribution maps of a single feature for various methods. On the left column, the input values of the satellite band B08 on March of the second year are visualized on the first map, along with the reference and predicted yield maps of a field from ARG-S dataset. The remaining columns contain the attribution maps of the same feature and field when different methods are used. The blue and red values distinguish between positive and negative attributions, which have contributed to either increasing or decreasing the predicted yield value for each pixel, respectively. We observe that only SVS produces a smooth attribution map, in contrast to the other attribution techniques where the results show pronounced spatial fluctuations. We verified that this pattern is consistent across the remaining features. Such behavior, indicative of high sensitivity, is undesirable because adjacent pixels often share similar input features and target values, implying that their attributions should also exhibit similarity [46, 47].

We further evaluate the choice of the baseline. While the results we described above are similar when using the padded baseline vector (see Appendix D), examining the attribution maps highlights a fundamental difference between the two options. In Figure 6, we select a field from the ARG-S dataset and visualize its attribution map for four satellite bands over three consecutive months. In Figure 6.a, the results using the mean baseline effectively differentiate between high and low yield regions by displaying positive and negative attribution values, matching the variance on the target yield map (refer to Figure 5). However, in Figure 6.b, the results of the padded baseline fail to capture similar patterns, with the maps consistently attributing values as either strictly positive or negative.

We would like to finally highlight an important advantage of using this method, which stems from game theory; Shapley values inherently accounts for feature interactions by approximating the impact of a feature when added to all possible subsets of the remaining features [48].

Based on our analysis of the evaluation results and the above discussion, we conclude that SVS is a reliable interpretation method for feature importance estimation, and we therefore conduct the remaining experiments with the SVS attribution method, using the mean baseline.

In this section, we investigate the interpretation results of the models trained on the different datasets and closely examine the spectral and temporal importance of each model. For each dataset, we examine the evaluation at each cross validation fold, and explain the model trained on the best-performing fold. We summarize the model evaluation results in Table IX in Appendix B, with the selected folds highlighted in green.

In the following, we connect the spectral importance results and the model performance when trained on vegetation indices. The normalized difference vegetation index (NDVI) is frequently used for yield prediction as highlighted in [12], and relies on the near-infrared (B08) and red (B04) bands. In the spectral importance results in Figure 7, we notice that B08 emerges as significantly important in all datasets, aligning with its universal relevance, while B04 is assigned high attributions in soybean datasets in particular. Additionally, the Sentinel-2 satellite imagery incorporates three red edge bands, of which only two are actually employed in the vegetation indices we utilized, as can be verified in Table X. However, the results underline that the third red edge band (B07) also exerts a high influence on the models across all experiments. Moreover, the first SWIR band (B11), also absent from the indices we used, has particularly high attributions in rapeseed and wheat datasets. This can potentially explain the substantial drop in the VI experiments for these two datasets, as reported in Table II and III.

The spectral importance enables the identification of the most important bands for each crop and region, and these results can efficiently be exploited to define specialized indices for the crop yield prediction task. Combining these results with the growth stage attributions, the vegetation indices can further be exploited to estimate the length of each growth stage.

In this section, we aim to compare our temporal importance results against well-established facts in crop growth and development research.

In soybean datasets, emergence, unifoliolate and trifoliolate stages collectively constitute the vegetation phase. The main stem nodes and their branching develop during this phase, influencing the canopy structure and the final number of nodes, as elucidated by Kumudini [54]. Since various environmental and genetic factors can impact the vegetative development, a poor canopy expansion can be identified at the last stage (i.e. 1st-25th trifoliolate) and linked by the model to lower yield values. The beginning pod stage is also strongly related to yield, with favorable temperature and moisture conditions resulting in a higher pod number per plant, more beans per pod, and larger seeds [43]. The following three stages, which extends from full pod to full seed stages, are important as the pods grow rapidly, and seed development begins, making it a crucial period for seed yield [43]. This seed-filling phase is proven to have a positive correlation with yield [55, 56], and any stress during this phase has a more significant impact on yield reduction than at other times[43]. Finally, the pronounced significance of the full maturity stage in Argentina can be attributed to the enhanced visibility of the pods. At this stage, approximately 95% of the pods on the main stem attain their mature pod color, offering a valuable visual indicator of yield, as supported by Kumidini [54].

One of the factors that could account for the differences in the results between Argentina and Uruguay is the notable influence of the flowering stage. As documented in [43, 54], soybean flowering is influenced by factors such as crop variety, day length, and temperature. Consequently, the variations observed in these regions may be attributed to differences in agricultural practices, including the application of fertilizers and herbicides.

In rapeseed, the attribution results of juvenile growth (i.e. from leaf development to heading) is attributed to the close connection between the seed yield and the number of pods per plant, which is supported by the dry matter produced during this period. There is a linear relationship between the cumulative production of dry matter until flowering and pod density, as documented by Diepenbrock [57]. He further highlights that the flowering stage is the most pivotal phase impacting rapeseed yield, which aligns with the high attribution of this stage in the bottom left plot of Figure 10. The ability of satellite data to capture this crucial stage is attributable to the increasing flower cover, which significantly enhances photon reflectivity and absorption to 60-65% of incoming radiation. At the fruit development stage, pod filling is initiated and continues through the ripening stage. Seed yield is linearly related to photosynthetically active radiation that is intercepted during this phase, which aligns with our attribution results. Moreover, while several factors during flowering can limit the yield, rapeseed has the potential for growth after flowering, which compensates for losses of buds, flowers and pods [57].

In the wheat dataset, significant importance is observed during the tillering stage, which is recognized to have great agronomic importance in cereals, as recognized by Acevedo et al. [58]. While it is the most important process governing canopy formation, it is also known that not all tillers produce spikes, which are the grain-bearing tip located at the top of cereal stems, as noted by Gallagher and Biscoe [59]. As shoot development progresses, spikes grow, the pseudo-stem becomes more erect, and leaf sheaths elongate. This results in a more distinguishable canopy density at satellite resolution, aligning with the high attribution of the shoot development stage, as illustrated in the bottom right plot of Figure 10. Additionally, florets are initiated during this stage; however, it’s noteworthy that less than half of these florets will complete anthesis [58]. As the crop progresses to the heading stage, the florets become ready for pollination and fertilization. It is at the flowering stage that the model leverages the blossoms for its predictions, as the number of flowers successfully pollinated directly influences the number of kernels per head. This observation accounts for the pronounced importance attributed to the flowering stage. During subsequent stages, the crop undergoes grain filling and moisture reduction during ripening, achieving maximum dry weight, and physiological maturity. During this period, the green color progressively fades, a change that can be captured by satellite data.

The analysis conducted above validates our temporal importance results aggregated over growth stages, strongly supporting that our models learned agronomically meaningful cues for accurate yield predictions.