Canopy temperatures distribution over soybean crop fields using satellite data in the Amazon biome frontier

ABSTRACT During the studied time window, between 2003 and 2010, there was an important increase of land use conversion into new soybean areas (first-time-use) in Mato Grosso state (MT) in Brazil. Uncertainties of future scenario of Brazilian agriculture and increase in the frequency of extreme events, such as the occurrence of high temperatures, is highly likely to produce yield loss on summer crops. The MT is the largest producer of soybeans and accounted for 28.2% of the national production in 2013. The objective of this study was to investigated specific characterization of land surface temperature distribution over the soybean crop fields canopies (canopy-LST) due to massive land use conversion into new soybean areas and its impacts on yield. Satellite imagery data from Aqua and Terra/MODIS sensors (Moderate Resolution Imaging Spectroradiometer) were compared to official agricultural statistics covering eight densely cultivated regions in the studied period. Results show that within the period from flowering to grain filling canopy-LST exhibits a non-negligible relation to yield. It is expected an additional loss of 4.9% on soybean yield for each 1oC of canopy-LST above the obtained optimal level of canopy-LST with 28.4oC, associated to the higher yield averages. The difference between overall average of canopy-LST and air temperature was found 4.2 oC.


Introduction
In the last two decades, the economic and environmental implications of increasing land area for soybean cultivation from deforestation and massive conversion of land use have become a major issue at the international policy agenda level (FAO, 2010;FAO 2011). Depending on future scenario of Brazilian agriculture which remains largely undetermined (Lapola et al., 2013;Cohn, VanWey, Spera, & Mustard, 2016), extreme meteorological fluctuations are able to cause severe effects on agricultural yields (Batistti & Naylor, 2009), causing an increase in environmental costs of food production (Godfray et al., 2010). Such changes can inhibit and adversely affect the crop development, since canopy temperature fluctuations often exceed the optimum range (Lobell & Burke, 2008;Gusso, 2013). Although seasonal meteorological fluctuations, such as high temperature and water stress occurrences, have been not always problematic, heat stress is increasingly gaining interest in scientific publications (Gusso, 2013), due to the need for assessment of future international demands (Lapola et al., 2013;Anderson et al., 2016) and market fluctuations (Batistti & Naylor, 2009) but especially for soybeans and maize with potential worldwide economic impacts by 2050 (UNEP -United Nations Environment Programme. Environmental Food Crisis, 2009. The Mato Grosso state (MT hereafter) is the largest producer of soybeans in Brazil and accounted for approximately 28.2% of the national production in 2013 (IBGE -Instituto Brasileiro de Geografia e Estatística, 2013). Although the total soybean crop areas did not increase drastically between 2003 and 2010, together with a decreasing trend of deforestation (Macedo et al., 2012;Lapola et al., 2014), study from Gusso, Ducati, and Bortolotto (2017b) has shown that there has been a massive conversion of other land uses into new soybean crop fields (first-time-use). Nowadays, soybean cultivation in Mato Grosso covers almost 8.0 million ha and the production level was more than 23.5 million metric tons in 2013 reaching a yield average of 3220 kg ha -1 , with 2930 kg ha -1 between 2003 and 2010 (IBGE -Instituto Brasileiro de Geografia e Estatística, 2013). From these expectations emerges the need for a better understanding of the impacts caused from recent environmental conditions on crop yield (Gusso, Arvor, & Ducati, 2017a). The objective of this study was to investigated specific characterization of land surface temperature distribution over soybean crop fields due to field homogeneity of massive land use conversion into soybean areas and its impacts on yield in the Mato Grosso (MT) state, Brazil.
Typically, EOS-MODIS (Earth Observing System-Moderate Resolution Imaging Spectroradiometer) satellite imagery data from NASA (National Atmospheric and Space Administration) have been applied to monitoring and modelling bioclimatic processes, crop cycle development and agricultural production in Brazil (Arvor, Meirelles, Dubreuil, Béegué, & Shimabukuro, 2012;Jasinsky et al., 2005). Although it is known that the grain number is especially affected by elevated temperature of vegetation canopy (hereafter canopy-LST) during flowering period (Batistti & Naylor, 2009), the quantitative assessment of production losses by measuring temperature impacts on yields is not well understood (Lobell & Burke, 2008;Schlenker & Roberts, 2009). Primarily based on the hypothesis that vegetation greenness and vigour levels are inversely related to canopy-LST, this study investigates recent spatially integrated relationship between canopy-LST from Aqua and Terra MODIS and official agricultural statistics of soybean yields.

Study area
Canopy-LST and yield data were analysed over densely occupied soybean crop areas in the MT state, between 2003 and 2010. The MT state (7.23°-17.87°S and 50.57°-61.52°W) soybean crop area covered almost 6.3 million ha in the crop year of 2010 (IBGE -Instituto Brasileiro de Geografia e Estatística, 2013). Cultivation in the MT state is mostly rainfed and the climate is tropical super humid (Af), with dry periods during the winter season (Köppen, 1948).
The sowing period for soybeans lasts from mid-September to early December (Arvor et al., 2012;Gusso et al., 2017a). Farmers plant soybeans after the onset of the rainy season, usually in October and soybeans remain the main crop, while maize or cotton is planted after the soybean harvest (Arvor, Jonathan, Meirelles, Dubreuil, & Durieux, 2011;Arvor et al., 2014). However, the duration of the rainy season is mainly related to the onset date (Arvor et al., 2014). The average annual precipitation is 1,610 mm with a standard deviation (SD) of 264 mm (INMET, 2009) but usually rain levels range between 1,200 and 2,000 mm (Arvor et al., 2011). Annual average air temperature is 24.4°C with SD of 1.5°C. During the period that covers phenological stage of flowering/grain filling, from December to February, the hottest month is January with 25.4°C (INMET -Instituto Nacional de Meteorologia, 2009).
The prevailing management practice is zero tillage farming (Plantio Direto in Portuguese), which is a low tillage planting and sowing practice greatly reducing soil erosion and organic matter degradation. Figure 1 shows the study area and the regional study limits established by IBGE (Instituto Brasileiro de Geografia e Estatística).

Data description and resolution
Data from several sources were used in this study as follows: (i) annual soybean agricultural statistics at the state and municipality levels (IBGE, 2013) were used to generate spatially interpolated maps of soybean yields on 2000 m spatial resolution and compare those to the state-scale results; (ii) annual soybean agricultural statistics at the state scale (CONAB, 2012) were used to compare and evaluate the results obtained from applied MODIS crop detection algorithm (MCDA) procedure for soybean area estimation and yield analysis; (iii) monthly accumulated precipitation data from the continually available meteorological stations distributed across the MT state (5) were obtained for the soybean crop period between September and December for each soybean crop year; (iv) enhanced vegetation index (EVI), MOD13Q1 product (collection 6morning overpassing) with 250 m of spatial resolution, in 16-day composite satellite images from MODIS/Terra data. MOD13Q1 product was used for DOY (day of year) 001 and 017 covering the entire MT with six tiles (H11V09, H11V10, H12V09, H12V10, H13V09, and H13V10). The MOD13Q1 product is initially required for crop profile identification and further crop area classification. Afterwards, it is also used for identification of vegetation canopy profile in order to generate the canopy-LST during maximum crop development; (v) LST, MYD11A2 product (collection 6afternoon overpassing) with 1000 m spatial resolution, in 8-day composite from MODIS/Aqua data. MYD11A2 product covering the entire MT using six tiles (H11V09, H11V10, H12V09, H12V10, H13V09, and H13V10) was used for the study period. Higher canopy-LST generally occurs during afternoon time. MYD11A2 products were combined two-by-two to obtain 16-day composites of maximum LST which were further combined with EVI at maximum vegetation development. Afterwards, both LST and EVI were retrieved from inside soybean crop areas in a way to be representative of the crop canopy-LST during the maximum vegetation development. The MOD13Q1 and the MYD11A2 products were combined for each one of soybean crop seasons from 2003 to 2010. The MODIS sensor provides an adequate imaging configuration for crop monitoring because of its almost-daily revisit rate combined with a spatial resolution of 250 m (Justice et al., 2002), which can be considered adequate for mapping large-scale agricultural fields in Brazil (Gusso et al., 2017a). Moreover, it has a good geometric quality that allows for time series and crop development analysis (Justice et al., 2002). The MODIS image data and products are freely available at https://lpdaac.usgs.gov/data_access/data_ pool; and (vi) administrative microregions from IBGE were selected based on the most intense production. Microregions are related to overall infrastructure characteristics as transportation and grain storage that it is supposed to have influence on crop management practices.

Crop area estimation
The use of a fixed calendar date to study remotesensing-based agricultural prediction models is not optimal (Bolton & Mark, 2013) because there is substantial variability in planting and maturing of soybean crops within the same general area (Brown, Kastens, Coutinho, Victoria, & Bishop, 2013). This is due to most recent crop management practices in the region as crop rotation, crop intensification and double crop system (Arvor et al., 2012;Spera et al., 2016). A double-crop system is used for soybeans in around 30% of agricultural areas in MT (Arvor et al., 2012).
In this study, soybean crop areas were identified by applying the MCDA which was developed by Gusso, Adami, Formaggio, Rizzi, and Rudorff (2012) for crop area estimation and tested by Gusso et al. (2017a) for the Mato Grosso ecoregions and prevailing conditions. MCDA is based on the identification of EVI values that are associated with soybean crop areas within period between late August and late October (DOY 225-289). By applying MCDA, it is possible to adequately mask canopy-LST image, soybean yield and spatially interpolated layers generated. Finally, it is possible to retrieve all pixels that are only associated with soybean crop area classification. At the municipality level, the MCDA estimates at 250 m resolution correlated well with the IBGE data, with an R 2 equal to 0.97 and root mean square deviation (RMSD) equal to 13,142 ha.
A crop area map was generated from a map composition that combined each crop year from 2003 to 2010, into one layer. The resulting crop area map covering 8.51 million hectares tagged all soybean crop pixels at a frequency equal to or greater than two events.

Correlation between vegetation index and yield
In MT soybean is usually planted from late September to early December, after onset of the rainy season, while the harvest period lasts from late January to early April (Arvor et al., 2012). Nowadays, most of local soybean cultivars have the total growth cycle ranging from 90 to 150 days (Sediyama, 2009). Using remote-sensing technology with MOD13Q1 data, Arvor et al. (2012) have found that 90-to 130day cultivars prevail in MT.
Studies from Gusso et al. (2012) and Gusso et al. (2017a) for monitoring the crop profile of soybean in Rio Grande do Sul and Mato Grosso found that vegetation indices correlate well with soybean yields. In MT, Gusso et al. (2017a) investigated the association of the prevailing physically driven conditions during DOY 001-033 to yield and the coefficients of determination ranged from 0.91 to 0.98, with overall result of R 2 = 0.96 (p ≤ 0.01). The period from DOY 001 to 033 (between January 1 and February 2) corresponds to the end of flowering period and the early grain filling stage for most areas in MT which is mostly sensitive to heat stress. In this way, by following the EVI profile, it is possible to identify the canopy-LST associated with the maximum EVI value, inside the time window of the flowering/grain filling period.
Yield maps for each crop year were obtained from the spatial interpolation 141 municipalities by using an ordinary kriging approach with spherical semivariogram, as shown in Figure 2(a). The kriging interpolation method is based on a variability trend in the distribution of sample values and the distance between them (Jian-Guo and Manson, 2009). Spatially interpolated maps of soybean yield data, obtained from official estimates of IBGE, were generated by allocating the municipality averages at the municipality centroid. Table 1 shows details of the spatial interpolation parameters used for each crop year.

Correlation between LST and yield
For earth science studies, the external manifestation of an object's energy can be detected by remotely sensed devices into an orbital path way using thermal scanning technology (Elachi, 1987). The retrieved LST obtained from on orbit sensors is a measure of surface temperature (also known as skin temperature) rather than air temperature (Elachi, 1987). The theoretical basis is Planck's Law of Radiation, which postulates that radiating energy from a black body, as predicted by Stefan-Boltzmann's Law, is distributed in the electromagnetic spectrum as a function of its temperature (Hecht, 1998).
When energy and water conditions are non-limiting for plant growth, considering complete canopy conditions, crops can present elevated canopy-LST due to aerodynamic resistances that suppress sensible heat transfer (Nemani & Running, 1997). In this way, considering that physiological activities of leaves are closely related to their actual temperature (canopy temperature) rather than to air temperature, LST can be used as a reliable measure of physiological activities of a vegetation canopy (Diak & Whipple, 1993). Furthermore, vegetative structures of soybeans which have a planophile canopy structure preserve spectral emissivity information as an approximation of a black body (Hecht, 1998;Salisbury, 1986) which is one of most important surface physical properties for accurate LST estimates. Under these theoretical circumstances, elevated canopy-LST can directly affect plant photosynthesis.
Among the major impacts of heat stress, its effect on photosynthesis is most important for yield due to inhibition of the crop growth rate (Schlenker & Roberts, 2009). Recent studies using satellite data series have also indicated that crop vulnerability to heat is due to plant damage and crop growth inhibition (Gusso, Ducati, Veronez, Sommer, & Silveira  Junior, 2014; Schlenker & Roberts, 2009), even when the average temperature reaches only one or two degrees above the optimal temperature for the culture (Sediyama, 2009). These effects are induced by reductions in canopy photosynthesis at the time of heat stress (Carmo-Silva et al., 2012), due to evaporative demand increases (Nemani & Running, 1997) or changes in the duration of the flowering period (Rodrigues, Didonet, Lhamby, Bertagnolli, & Luz, 2001). Optimal canopy-LST values between 298 K (25°C) and 309 K (35°C), with values nearby 303 K (30°C) are most frequently considered the ideal conditions for soybean development (Schlenker & Roberts, 2009). Through a different ecosystem than MT, an important contribution regarding this research is given by Gusso et al. (2014) when studying the southern region of Brazil. Those studies observed that slightly warmer than the optimal conditions for growth can potentially increase drought effects (if any) and yield loss, by overheating the vegetation canopy and intensifying plant damage. Evapotranspiration reduction leads to more energy available for sensible heating of the surface and canopy-LST ultimately, which is partly controlled by soil moisture availability (Sandholt, Rasmunssen, & Andersen, 2002). As predicted by literature, they observed that heat that occurs during flowering and pod formation (R1-R5) affects grain number, which is closely related to the yield (Gibson & Mullen, 1996) instead stress occurred during grain filling (R5-R7) which reduces the grain size (Gibson & Mullen, 1996). The authors obtained in the southern region a correlation R 2 equal to 0.82 with RMSD = 14.4% between canopy-LST and yield at microregional scale during February. Thus, the following hypotheses have been proposed: (1) different spatial configurations of canopy-LST from one region to another leads to different levels of regional production; (2) vegetation greenness and vigour levels are inversely related to canopy-LST; (3) it is possible to draw an association between the canopy-LST and the soybean yield by means a continuous mathematical function; (4) the period to which soybean vegetation reaches the maximum canopy coverage (MCC) can be associated with phenological stage of flowering/grain filling; (5) vegetation index profile, between January 1 and February 2 can be used for locate the MCC and so, the canopy-LST associated with this moment.
It is important to note that between 1991 and 2010 it is observed a rapid increase of yield. In this period, state yield average raised from 2,370 to 3,015 kg ha -1 , with an increasing rate of 38.6 kg ha -1 year -1 , due to technological and logistics improvements as well as crop breeding and agricultural management practices. Even so, because this study evaluated the spatially distributed characteristics of canopy-LST averages from 2003 to 2010, it did not perform any adjustment for yield due to technological improvements through time.

Classification procedure
In order to match the LST and EVI products in a way that one can obtain LST as representative of vegetation canopy during the period of MCC, imagery data were rearranged and grouped them according to the following procedure: (1) the MOD13Q1 product already represents a 16-day composite of maximum EVI. Further calculations were performed within 16 days of EVI, starting from DOY 001 which is associated with the flowering period for soybean. (2) To retrieve the most representative LST data for MCC, the first step was to obtain the maximum LST composition over 16 days by combining two MYD11A2 products (an average of 8 days) that cover the same 16 days DOY period of MOD13Q1 in item 1. (3) After that, steps 1 and 2 were also performed for the subsequent 16-day composites (DOY 017) of MYD11A2 and MOD13Q. By doing this, it is obtained two subsequent 16 days composites of maximum LST and two subsequent 16 days composition of maximum EVI. (4) By combining the two LST composites of 16 days (DOY 001-017) as a function of the maximum EVI, we obtain a 32-day composite of canopy-LST when crop canopy vegetation is its maximum. This is an important step because the high incidence of cloud cover exerts a confusing factor during key periods of the vegetation profile (DOY 353 to 033) in MT (Gusso et al., 2014(Gusso et al., , 2017. Additionally, in order to overcome the cloud cover and cloud shadows challenge, the 8-day product (MYD11A2) was combined into a maximum LST of 16 days. The resulting image is a 32-day composite that represents the crop canopy-LST that occurred when EVI was at maximum between DOY 001 and 032. Analysis of earlier periods of canopy-LST (prior to DOY 001) was not performed due to the heterogeneous surface characteristics which predominate during the initial phenological stage. Such heterogeneous conditions in the surface can induce inaccuracies of several degrees in the near-surface air temperature and consequently, in the canopy-LST as well, because of the absence of the soybean vegetation canopy in the initial stages of crop profile development. (5) Distribution of canopy-LST imagery was generated for the phenological stage covering a time window of 32 days, starting from DOY 001 (as presented in Figure 2a). (6) The resulting three levels of information, spatially interpolated yield maps, canopy-LST and EVI imagery, were recombined with a spatial resolution of 250 m. (7) Finally, it is obtained the averages of the most intensive agricultural microregions which are defined by the weighted crop area and the eight microregions from IBGE. As defined at section Crop Area Estimation, from the total 8.51 million hectares, only the selected eight densely occupied regions of soybean crops were included in this study, which covers 6.00 million ha in MT, corresponding to 70.5% of total soybean area.

Satellite measurements
After obtaining subsets of soybean crop areas, all data were compared and analysed at the state and regional subset scale, as shown in Figure 2. The groups were selected based on the regions to generate state subset averages of canopy-LST for each crop year. Figure 2 (a) shows the averages maps of canopy-LST and Figure 2(b) shows spatially interpolated map of soybean yields between 2003 and 2010.
The LST data were analysed on the Kelvin scale, and the yield is expressed in kg ha −1 . The spatial yield variability was analysed by combining it with canopy-LST and a vegetation index at the regional scale, based on the assumption that the agrometeorological conditions are evenly distributed across the entire state. The distribution of canopy-LST showed occurrences of lower temperatures in the central-west region (dark red in Figure 2(a)) where soybean crop fields with higher regional yield averages prevailed (light blue in Figure 2(b)). Accordingly, the distribution of canopy-LST with higher temperature occurrences was observed in the south-eastern region (light blue in Figure 2(a)) where soybean crop fields with lower regional yield averages prevailed (yellow in Figure 2(b)).

Data retrieval
The first-order polynomial regression analysis was performed for the averages of regional subsets obtained by comparing canopy-LST (independent variable) and spatially interpolated yield maps (dependent variable) for each crop year between 2003 and 2010. Better correlations were found by using a first-order polynomial function between canopy-LST and the yield under the general mathematical formula: where y is the dependent variable given by yield; β is the slope which represents the increment of yield in response to changes in canopy-LST; α is the intercept. When β is negative, α can be projected as the minimum "x" to obtain "y = 0".
Microregional scale: spatial analysis At microregional analysis, absolute yield averages for all selected regions showed an inversely related trend (Figure 3) between yield and canopy-LST, as predicted by the theory, with R 2 equal to 0.44 and RMSD of 153.2 Kg ha −1 (p ≤ 0.05). Overall averages at microregional scale, from recent imagery analysis (Table 2) are canopy-LST 302.8 K (29.6°C) with SD of 2.62 K and yield 2702.5 Kg ha −1 with SD of 353.9 Kg ha −1 . The difference between air temperature in the studied period (25.4°C) and overall average of canopy-LST (29.6°C) was found 4.2°C.
The higher yield value was obtained in Alto Teles with 3020.1 kg ha −1 and canopy-LST 301.6 K (28.4°C) which is the optimal level. The optimal level can be interpreted as the canopy-LST average associated with the higher yield averages. The lower yield averages occurred in Rondonópolis with yield 2341.3 kg ha −1 and canopy-LST 304.0 K (30.85°C). Table 2 summarizes overall results obtained from the values presented in Figure 3.
Regarding spatially distributed extremes, Table 2 shows the minimum and maximum extremes at pixel scale. Minimum yield occurrence was in the microregion of Rosário do Oeste with 1092.5 kg ha −1 when SD was 542.5 Kg ha −1 and the maximum yield Minimum extreme canopy-LST occurrence was in the microregion of Parecis with 295.7 K (22.5°C) when SD was 2.33 K and the maximum extreme canopy-LST occurrence was in Canarana with 311.7 K (38.55°C) and SD was 2.47 K.
Microregional scale: temporal analysis Considering the regional subset scale at yearly averages the parameters β and α were evaluated according Equation (1). As evident in Figure 4, the regression analysis shows that all the parameters β have a negative declination. Table 3  Although the canopy-LST effects through time has not been extensively studied here it is possible to see either from Table 3 and Figure 4 that in the period from 2003 to 2010 the better yield occurrences are for almost all years, closely associated with the microregion of Alto Teles. There is only one year (2003) which the optimal yield also occurred with a higher canopy-LST average.
As expected from theory, the well-known negative mathematical relationship between vegetation index and canopy-LST on local and regional scales is also observed for yield and canopy-LST. Considering the size of each of microregional area studied and the number of pixels associated with 2000 m, the obtained differences from Table 2 and Figure 3 between optimal and poorest yield levels (648 Kg ha −1 ) and the related canopy-LST (2.4 K) are consistent.
Overall results obtained in this study, as summarized in Table 2 and Figure 3, indicate that there is an evident condition of lower yields associated with higher canopy-LST values. This means that the spatially distributed yield in DOY 001 (between January 1 and February 2), which corresponded to the end of the flowering period and the early grain filling stage, was inversely related to canopy-LST. Even when analysis was performed on yearly basis there has been an overall microregional similarity of the linear regression trends. In this way, the closely similarity between Figures 3 and 4 is due to the existing prevailing microregional characteristics that relates canopy-LST and yield averages. This clearly point out that there is a spatial dependence of yield to prevailing physically driven conditions as temperature.
As shown previously in Table 2, the microregion of Alto Teles has exhibited the optimal canopy-LST conditions while Rondonópolis the non-optimal canopy-LST conditions. Study from Spera et al. () also evaluated impacts of different factors as temperature parameters, related to expansion and abandonment of agriculture areas, and found that agriculture in MT may eventually become concentrated on the most favourable areas where prevails better conditions. However, they use climatic air temperature database from Climate Research Unit (V.3.21) temperature data at 50 Km resolution for 1980-2000 merely based on statistical analysis.
On the microregional subset scale, higher soybean yield averages are also associated with lower canopy-LST occurrences in the MT, which also agrees with the results obtained by Gusso et al. (2014) for crop fields in Rio Grande do Sul ecoregions, Brazil. Thus, even general results from imagery analysis point out that the influence of canopy-LST on yield is nonnegligible.
Considering the equation exhibited in Figure 3, if canopy-LST tends to occur slightly above the actual averages, obtained from Table 2 (29.6°C), with roughly 2°C, there will be a constant loss of 9.9% on soybean yield, decreasing from an average of 2702 to 2441 Kgha −1 .
Linear regression approaches from the relation of LST and vegetation index (Nemani et al., 1993;Sandholt et al., 2002) show that higher declinations with steeper linear regressions are representative of heat-induced water stress conditions among selected microregions. This is because higher constant β values can be seen as trend for either spread of canopy-LST and yield distribution. On the other

Conclusion
This study evaluated the spatial distribution of canopy-LST and soybean yields on regional scale. Although recent technology development and management practices have increased the potential for higher soybean yield averages, lower canopy-LST values during flowering/grain filling still result in better yield averages in MT. Although this study did not explore specific agrometeorological causes under which canopy-LST affects soybean yield and production, this study clarifies that the characterization of canopy-LST distribution over soybean crop fields are associated with a non-negligible influence of its spatiality in MT. Future research must investigate the relation between specific crop development conditions to the observed agrometeorological details.

Disclosure statement
No potential conflict of interest was reported by the author.