Estimating and understanding crop yields with explainable deep learning in the Indian Wheat Belt

The performances measured as NSE for different models and input datasets are compared for all years averaged and for one year separately (2012) in table <a xmlns:xlink="http://www.w3.org/1999/xlink" class="tabref" href="#erlab68act2"="">2</a>. The results for 2012 were shown to demonstrate that the overall good performance among all years does not necessarily guarantee accurate predictions for abnormal years. 2012 was a peculiar year characterized by very high yields compared to average (figure S8), as indicated by poor performance of the null model (table <a xmlns:xlink="http://www.w3.org/1999/xlink" class="tabref" href="#erlab68act2"="">2</a>). CNN models provided the best results overall (best NSE of 0.868 among all years in the 'Incl. district' setting, compared to the best NSE of 0.757 for RR and 0.836 for RF) that were stable among model runs (S4), which shows applicability of such models to dense time series for yield estimation. The null model was already a good predictor for the yield among all years (NSE of 0.812) and it performed better than any model that excluded district information, as the majority of the yield variation came from the spatial variation (S5). However, when considering the abnormal year 2012, the null model performed much worse than the other models (NSE of 0.288), which demonstrates the importance of using satellite observations to estimate yields for atypical years.

Performance of CNNs improved with the increase in the number of input parameters, both in terms of the best performing models, as well as among all tested model settings (table <a xmlns:xlink="http://www.w3.org/1999/xlink" class="tabref" href="#erlab68act2"="">2</a>). When testing models with one or two variables to evaluate the parameters importance, the model performance and variable ranking were not consistent among the years, with 2012 showing the worst performance (figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf3"="">3</a> and <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf4"="">4</a>). The models performed much better and had a smaller performance variation if the information on the district was included. However, in 2012, some model combinations gave better results without district information, which is related to the fact that the null model performed poorly for this year.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In case of models with one variable only and no district information ('No district' in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf3"="">3</a>), VIs performed better among all the parameters (and NDWI the best). Models with day-length as the only variable had a very good performance, even though this parameter does not carry any information about the crop condition. It suggests that the model was actually trying to learn the specific day-length patterns related to districts, as yields have a clear spatial pattern in this region (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf1"="">1</a>) that is somewhat similar to the latitudinal distribution of day-length (figure S9).

When the district information was fed into the model ('Incl. district' in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf3"="">3</a>), the best performing model across all years was the one using NDWI, followed closely by SW<sub="">down</sub> and other VIs. The good performance of NDWI and SW<sub="">down</sub> was reproduced for models with two variables (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf4"="">4</a>). However, this behavior was not consistent among all the years for neither one- nor two-variable models. Good performance of NDWI and SW<sub="">down</sub> in 2012 (figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf3"="">3</a> and <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf4"="">4</a>) could be a leading factor responsible for their high overall ranking, as the differences among the models in this year were much larger than for other years. For example, in case of 2002, the ranking of pairs including NDWI and SW<sub="">down</sub> were among the worst, but the general variability was very small. These results suggest that accounting for various parameters is important, as the conditions that limit or boost the crop yields vary among the years.

Overall, the main features of RAMs were the same for all the ensemble members (figure 5), as well as for the models re-trained on all the input data (CNN${}_{2,\mathrm{all}}$ and CNN${}_{10,\mathrm{all}}$). We compare in detail RAMs of the district Patiala for 2006, for selected CNN models using two (CNN${}_{2,\mathrm{all}}$) or ten (CNN${}_{10,\mathrm{all}}$) input variables in figures 6(a)–(b). As mentioned before, the average of RAM plus bias (bias both of the whole model and the district-specific, see S6), equals the estimated yield. Even though RAMs are products of the complex nonlinear interactions of the input data, some basic patterns can be directly inferred. For example, the overall shape of RAMs is closely related to the growing cycle as shown by NDWI. In case of the 2-variable model (figure 6(a)), all observed temporal patterns of RAM can be related to NDWI or T_min, as these were the only input parameters. As a result, a small dip in RAM around day 100 can be associated with the increase in T_min, as NDWI during this time was steadily increasing, which suggests the negative impact of higher T_min on the crop yield. To support this claim, we modified T_min and calculated the resulting changes in RAMs as shown in figure 6(c). Two small changes were performed on the T_min data—in the first case we removed the T_min peak around day 100 (shown in blue in figure 6(c)) and in the second case, we increased T_min around day 150 (shown in yellow in figure 6(c)). The changes in the resulting RAMs—which directly relate to changes in the estimated yield—are highlighted in blue and yellow in figure 6(c). The decrease in T_min led to the increase in RAM and removed the previously observed dip. On the other hand, increase in T_min led to a similar decrease in RAM, and therefore a decrease in the yield. However, analogous T_min modifications for CNN${}_{10,\mathrm{all}}$ resulted in different magnitudes of changes in RAMs due to decrease or increase in T_min. The increase in T_min towards the end of the growing season had a much bigger impact on RAM than the decrease in T_min earlier in the growing season, which led to a positive but very small change in RAM.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To better identify the parameters and features that led to certain responses in RAMs, we compared selected input variables and RAMs for Patiala for two very different years: the year 2006, when the yield was quite poor (4.233 t/ha), and 2012, when the yield was very good (5.473 t/ha) in figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf7"="">7</a>(a) and (b). In 2006, RAM demonstrated overall lower values and many strong drops, while in 2012, RAM showed consistently high values despite not much higher VIs. To analyze the drivers of this variability, we exchanged, one in turn, four input variables (<em="">T</em><sub="">min</sub>, <em="">T</em><sub="">max</sub>, SW<sub="">down</sub>, precipitation) among 2006 and 2012 to check how this would affect changes in both RAMs and estimated yields (figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf7"="">7</a>(c)–(j)). For example, SW<sub="">down</sub> was mostly responsible for creating the dips in RAMs in 2006 for Patiala—as many of them were filled when SW<sub="">down</sub> from 2012 was used, especially the big drop around day 160 (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf7"="">7</a>(g)). Overall, the yield in 2006 increased using any of the input parameters from 2012, but the highest increase was observed for <em="">T</em><sub="">max</sub>, followed by SW<sub="">down</sub>, <em="">T</em><sub="">min</sub> and only slightly for precipitation. The same rank of variables in terms of the magnitude of the impact was observed in the case of modifying the input variables in 2012.

Zoom In Zoom Out Reset image size

Similar features were consistent among all districts in 2006 and 2012, as shown by anomalies (positive in red, negative in blue) in RAM, NDWI, <em="">T</em><sub="">max</sub>, <em="">T</em><sub="">min</sub> and SW<sub="">down</sub> in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf8"="">8</a>. In general, these years were relatively sunny, but in 2012, both <em="">T</em><sub="">max</sub> and <em="">T</em><sub="">min</sub> had strong negative anomalies during the second half of the growing season. The rapid and short negative anomalies in RAMs were usually related to negative anomalies in SW<sub="">down</sub> (highlighted with green boxes), which sometimes were also accompanied by positive anomalies in <em="">T</em><sub="">min</sub> (as cloudy skies capture the long-wave radiation emitted by Earth at night). It is noted that this link may lead to a wrong interpretation of the effect of <em="">T</em><sub="">min</sub>, when considered alone (e.g. a model with <em="">T</em><sub="">min</sub> as the only weather input). Positive anomalies in RAM in 2012 were connected to strong negative anomalies in <em="">T</em><sub="">max</sub> and <em="">T</em><sub="">min</sub> (highlighted with magenta boxes).

Zoom In Zoom Out Reset image size

We examined if our CNN model could be potentially useful for predicting the yield during the growing season by analyzing the impact of shortening the input time series in 25-day steps on the model performance. As compared to the null model, CNN<sub="">10</sub> had better prediction from day 145, which corresponds to the last week of February (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf9"="">9</a>), In general, the month of February is the time of the anthesis stage, conditions during which are crucial for the final crop yield. Therefore, covering this time period is essential for a good yield estimation. From day 195 (middle of April), which widely corresponds to the harvest time in this region, the model performance did not increase. This suggests that after harvest the additional data were not useful anymore for the yield estimation. Although this may seem trivial, it is nevertheless important that model performance appears insensitive to post-harvest data, as adding such data could occur whilst in an operational context.

Zoom In Zoom Out Reset image size

Even though our DL modeling approach was constrained by the requirement for network explainability (global pooling layer after the convolutional layers to facilitate RAMs), it outperformed RF and RR in the yield estimation task, though RF provided a strong performance improvement as compared to RR. Although DL models are already widely used for a variety of problems in image and signal processing domain (LeCun <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib19" id="fnref-erlab68acbib19"="">2015</a>), DL applications for crop yield estimation are still rare (e.g. You <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib39" id="fnref-erlab68acbib39"="">2017</a>, Crane-Droesch <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib4" id="fnref-erlab68acbib4"="">2018</a>), as the transfer of available tools from the field of ML into the environmental applications requires adaptations that account for specificities of the data. Our results demonstrate that CNNs can be a valid tool for capturing complex processes in agricultural systems, and their application can be therefore further explored for other crops and regions.

To generalize well in the yield estimation problem, it is important to train models over several years and simultaneously account for multiple vegetation and environmental variables, as the main yield drivers vary among seasons. Since many of the input variables in the crop yield models are correlated and inter-related, multiple parameters should be considered when trying to understand the impact of variables on the crop yield in the model. For example, the impact of <em="">T</em><sub="">min</sub> varied depending on whether it was the only input meteorological variable or one of many (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#erlab68acf6"="">6</a>).

RAM's general shape reflected VIs, which emphasizes the importance of the length of the growing period and agrees with previous studies (Lobell <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib22" id="fnref-erlab68acbib22"="">2012</a>, Jain <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib17" id="fnref-erlab68acbib17"="">2017</a>). Meteorological variables used together with VIs are expected to have only a marginal contribution to yield estimation, as their impact is already reflected by vegetation growth that is captured by VIs. Overall, SW<sub="">down</sub> turned out to be the most important meteorological variable for yield estimation, even though the ranking varied significantly among the years. Comparison of performances among models with one or two input variables, as well as the RAM analysis, showed that high yields in 2012 were associated with low temperatures accompanied by sunny conditions during the growing period, though most recent studies focus primarily on the impact of increasing temperatures in the Indian Wheat Belt (e.g. Duncan <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib5" id="fnref-erlab68acbib5"="">2014</a>, Jain <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib17" id="fnref-erlab68acbib17"="">2017</a>, Song <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib33" id="fnref-erlab68acbib33"="">2018</a>) and in other regions (Lobell <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib20" id="fnref-erlab68acbib20"="">2005</a>). As an important driver of photosynthesis, SW<sub="">down</sub> affects the resources that crops can build during the growing season, which is apparently not fully reflected in the VIs. High radiation levels are also especially important during the critical period for grain number determination (20–30 days before anthesis to ten days afterwards) (Fischer <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib9" id="fnref-erlab68acbib9"="">1985</a>). Analysis of RAMs showed that the model was trained to recognize events of decreased light as having negative impact on the crop yield, which suggests that the decline of photosynthesis due to decreased radiation was larger than the benefit due to increased diffused light. This agrees with the negative effects associated with a reduction in sunlight shown on a global scale (Proctor <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib28" id="fnref-erlab68acbib28"="">2018</a>). Decreased penetration of sunlight for the crop has already spurred demand for breeding for high photosynthetic or radiation-use efficiency (Joshi <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib18" id="fnref-erlab68acbib18"="">2007</a>). These observations also suggest that sun-induced fluorescence could be a good direct proxy for the crop yield, as it carries information on light conditions, and was previously shown to perform better than VIs for yield estimation (Song <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib33" id="fnref-erlab68acbib33"="">2018</a>). Although NIRv was found superior to other VIs for estimating GPP (Badgley <em="">et al</em> <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#erlab68acbib2" id="fnref-erlab68acbib2"="">2017</a>), we did not detect any benefit of using it for the yield estimation.

The novelty of our approach extends beyond obtaining a satisfactory model performance, as we adapted the activation mapping to the regression problem for the temporal data in order to localize important patterns for the yield estimation. The analysis of RAMs applied in the time dimension facilitates DL network interpretability and provides the needed transparency on what DL models learn and how they accomplish their prediction. Analysis of the different patterns in RAMs can provide a consistent description of how the network captures the weather effects in relation to the temporal progression of the crop, and can ultimately lead to a valid physical interpretation. The natural way to interpret RAMs is to regard it as a kind of derivative of the yield function in time. In the global pooling layer this function is then integrated out by the network to estimate the crop yield, which as a result reflects the accumulated effects of crop growth and meteorological conditions. Thus, the application of RAM allows for an instinctive analysis of how certain events in the time domain impact the estimated crop yield and enables a comparison with known drivers. In our case, RAMs were roughly reflecting the development of the crop and at the same time were sensitive to the variability of meteorological conditions, which somehow reflects the basic approach of relating yield to the accumulated biomass. Such an approach can increase the confidence in the model (as it captures the processes that are expected to be relevant), but might also draw attention to so far neglected features or model weaknesses. For example, although high temperatures leading to low yields are often considered in the Indian Wheat Belt, the importance of sunny conditions for good yields is rather neglected. On the other hand, RAMs varying towards the end of the analyzed time period might suggest that the model has not completely learned to ignore the time after the harvest. In general, the proposed methodology facilitates the application of DL in agriculture, not only to improve yield estimation and prediction but also to gain insight into the key drivers of crop yield.