Atmospheric Effects on Spatial Crop Yield Modelling Using Landsat’s Imagery

Field Spectroscopy

Field spectroscopy is a technique used to measure the spectral characteristics of ground surfaces in the natural environment [1]. Field data at or near ground level permit the allocation of points or areas on satellite imagery to corresponding regions on the ground and therefore enables the retrieval of information on the spectral signature of different features [2]. In this study, a field spectroradiometer, GER 1500 was employed, with a range of spectrum from 350 nm to 1050 nm. The general and calibration characteristics of field spectroscopy described by Dekker AG, et al [3]. Spectral signatures can be used to discriminate between healthy and non-healthy vegetation [4, 5, 6, 7, 8].

The role of field spectroscopy has been discussed by Milton and McCoy KR [1, 9]. Field spectroscopy is used to characterize the spectral behavior of ground targets and monitor their suitability as appropriate targets over time. It helps in the inter-calibration of data from different platforms or the calibration of data from the same sensor at different times. Satellite images are used to infer reflectance with and without atmospheric corrections and compare it to the reflectance of a spectroradiometer (which is considered as corrected from atmospheric effects). Ideally, the application of an algorithm for atmospheric compensations, on an image, should finally provide the corrected reflectance which is the reflectance provided from the spectroradiometer.

Crop Yield Modelling

Quantifying crop production at regional scales is critical for a wide range of applications and remote sensing offers great potential for monitoring regional production [10]. Generally, knowledge of crop phenology is combined with multi-temporal imagery to estimate Durum wheat crop yield in the area of interest. Semi- empirical models have been developed to explore crop yield predictions for many different places, including the one used in this paper, which provide valuable information regarding the spatial and temporal distributions of crop production [11, 12]. Figure 1 shows the plots that were used to infer the necessary data.

Click to enlarge

Crop production is strongly related to many crop, meteorological and environmental factors. Field experiments from the early 1960s have shown that these factors proportionally affect the production of a crop while new experiments have indicated the possibility to estimate a crop’s yield by combining these factors along with remotely sensed data [11, 12, 13, 14]. Therefore, the purpose is to correlate basic factors that affect crop production of Durum Wheat using remote sensing and statistical analysis. The use of remote sensing has been proven to be effective in monitoring the growth of agricultural crops and in irrigation scheduling and efforts have been made to develop various indices for different crops of different regions throughout the globe. The production and prediction of crop yield have a direct impact on year-to-year national and international economies and play an important role in the food management [15].

The semi-empirical model used by Agricultural Research Institute (ARI) for estimating Crop Yield in Cyprus is described in Table 1. It has to be mentioned that the model is best fitted when applied in small scale analysis and when using medium to high resolution satellite imagery (5-30m) where is very effective as mentioned by [16].

$$ \begin{array}{l} \mathbf {x} 1 = N D V I \\ \mathbf {x} 2 = L A I \\ \end{array} $$

x3= Soil humidity (depth of 40cm)

x4= Crop temperature The model has its best results when applied at the stage of grain filling (stage duration: 8-12 days) with a very high R2 of 0.886 which is a high level of correlation indicating that the independent variables can explain by 88% the dependent variable with standard error of estimate of 0.12 which is very low. The F statistic value is 9.085 indicating that the overall regression model has a good fit for the data since the observed F (9.085) value is higher than the statistic (F4,45=5.7) for p<0.05 illustrating that the hypothesis (Yield can be predicted from NDVI, LAI, Soil Humidity and Crop Temperature) is valid. Significance of 0.002 (less than 0.05) suggest that, overall, the model can statistically significantly predict the outcome variable (Crop Yield). Thus, the prediction model is very accurate and crop yield can be calculated or estimated from the independent parameters NDVI, LAI, SH and CT.

Literature Review

The effects of the atmosphere on spectral signatures and vegetation indices have become an important issue in relevant scientific literature since the 1980s [17]. Atmosphere is a primary source of noise for accurate measurement of surface reflectance with optical remote sensing [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. Atmospheric effects are a result of molecular scattering and absorption, and influence the quality of the information extracted from remote sensing measurements, such as, vegetation indices. Such errors, caused by atmospheric effects, can increase the uncertainty up to 10%, depending on the spectral channel [29]. Hadjimitsis et al. highlighted the importance of considering atmospheric effects when several vegetation indices, such as NDVI were applied to Landsat TM/ETM+ images for agricultural applications [30]. In their study a mean difference of 18% for the NDVI was recorded before and after the application of darkest pixel method Figure 2.

Click to enlarge

Therefore, removal of the atmospheric effects is an important pre-processing step required in many remote sensing applications, since it is needed to convert the at- satellite spectral radiances of satellite imagery to them at- surface counterparts [31, 32]. Moreover, if an image is to be used for change detection and monitoring purposes, it is essential that adequate atmospheric and radiometric processing be applied, in order to bring all scenes to a common radiometric datum.

Several atmospheric correction algorithms are found in the literature, ranging from simple to sophisticated [33, 34, 35]. Hadjimitsis et al. classified these algorithms into the following two categories [35]:

• Category (A): Absolute corrections (corrections that lead to surface reflectance). This category can be subdivided into two sub-categories: image-based atmospheric corrections (for example, Darkest Pixel, Covariance Matrix method) and corrections using independent data for atmospheric optical conditions (including in situ measurements or historical records) using physical-based algorithms;
• Category (B): Relative corrections (corrections that do not produce values of surface reflectance).

The problem of atmospheric effects is especially significant when using multi-spectral satellite data for monitoring purposes such as agricultural or land use studies [5, 16]. Hence, it is essential to consider the effect of the atmosphere by applying a reliable and efficient atmospheric correction during pre-processing of digital data. Atmospheric corrections interfere at the stage of pre-processing after Digital Numbers (DN) is converted to Radiance. At this stage atmospheric correction algorithms are applied in order to remove the effects from the atmosphere. A considerable investigation on the effect of the atmosphere on dark targets, in the area of interest, has been already examined; the impact of atmospheric correction on vegetation indices as well on crops as separate targets has been raised and is investigated in this paper [5, 35, 36, 37, 38].

The advantage of the image-based algorithms is the fact that they do not require any ground data during satellite overpass [39]. One of the simplest, fully-imaged based, is the Darkest Pixel (DP) algorithm, categorized in Category (A). The DP principle is based on the assumption that most of the signal reaching a satellite sensor from a dark object, such as, deep inland water bodies (as in the case of the Asprokremmos Dam followed in this study) is contributed to by the atmosphere at visible wavelengths and on the assumption that the near-infrared and middle infrared image data are free of atmospheric scattering effects [39]. Therefore, the pixels from dark targets are indicators of the amount of upwelling path radiance in that band. The atmospheric path radiance adds to the surface radiance of the dark target, giving the target radiance at the sensor. The surface reflectance of the dark target, for the DP algorithm, is approximated to have zero surface reflectance. A modified adaptation of the DP method is to assume a known non-zero surface reflectance of the dark target based on ground truth data (e.g., spectroradiometric measurements).

Objectives

Main objective of the study is to check whereby atmospheric effects can act on yield prediction through remote sensing, a method which is walking over the classic direct methods. Many atmospheric correction methods have been proposed for use with multi-spectral satellite imagery [40]. Such methods consist of image- based methods, methods that use atmospheric modeling and finally those methods that use ground data during the satellite overpass. International literature [2, 7, 41] always premises the need for atmospheric correction on satellite images, especially when it comes to prediction models.

Hence the rhetorical question raised in this study, is “what could happen if atmospheric corrections are not taken into account or if an atmospheric correction algorithm is not working properly”. Never before scientists that apply modelling with the ‘spatial’ parameter have indicated that they use effective atmospheric removal algorithms. The usual path is to apply the model or algorithm onto images of indifferent spatial scales and platforms (for example satellite or UAV image). This paper comes to support modellers and introduce a sound novel methodology for more accurate results. To assist this procedure, authors have used a handheld field spectroradiometer to collect the reflectance of Durum wheat at the time of the satellite overpass in order to compare the correct reflectance (from spectroradiometer) and the reflectance of the satellite images without applying any atmospheric corrections. Based on this the authors have proceeded to estimate vegetation indices, LAI and Crop Yield with and without any atmospheric corrections. Finally, these pairs were compared to find out what could happen to them if atmospheric corrections are not applied to satellite images.

The Study Area

Durum wheat is a vital crop for Cyprus. The main wheat types are the Ourania and Ekavi which are the local types developed by the Agricultural Research Institute of Cyprus. Thus, prediction of wheat production can be useful to policy makers in order to establish policies regarding the market economics. In this study, 25 plots were cultivated with Durum wheat in the district of Paphos province, which were then marked and observed for this study. The area of interest is in the southwest of Cyprus, in a coastal strip between Kouklia and Yeroskipou villages, in the Akhelia area of Pafos. The area is a coastal plain with a seaward slope of about 2% and consists of deep fertile soils. The area is divided by three major rivers, the Ezousa, Xeropotamos and Diarizos. The study area is a traditionally agricultural area with annual and perennial cultivations and is irrigated by Asprokremnos Dam, one of the biggest dams of Cyprus.

The area is almost at sea level (altitude 15 m) and is characterized by mild climate which provides the opportunity for early production of leafy and annual crops. The uniform and moderate temperatures, attributed to the permanent sea breeze of the area, and the relative humidity, are conductive to the early production of fruits and vegetables. Cereals are also cultivated in the area. A typical Mediterranean climate prevails in the area of interest, with hot dry summers from June to September and cool winters from December to March, during which much of the annual rainfall occurs with an average record of 425mm. The whole are is generally homogenous in terms of soil and morphology while the dominant type of soil is the Cambisol (calcic and chromic types).

Wheat in the Akhelia area, Pafos, is sown in late November to mid-December using a sowing machine. In general, the study area is rainfed but, under severe drought conditions, is irrigated. The fields are treated with a basic 20:20:0 fertilizer upon sowing and again at a later time receive an additional top fertilization using a 34.5:0:0 formula. After the 4th leaf, the study area is treated with a wheat-triticale specific pesticide that controls both broad and narrow leaf weeds. The model application is taking place at the grain filling stage which is normally in the beginning of April and ends in the mid- April.

Overall Methodology

The overall methodology lies in two levels: a) estimating Vegetation Indices along with LAI and Crop Yield using satellite images, field spectroscopy techniques and the Crop Yield prediction model and other instruments and b) comparing using statistical techniques these factors (VI, LAI and Crop Yield) to the real values for inferring the effects of atmosphere on the procedure for predicting Crop Yield using remotely sensed data . More analytically, the methodology can be categorized in the following four steps (a) to (d): a) A field campaign, using a field spectroradiometer, was undertaken in the area of interest (25 plots cultivated with Durum wheat) during April 2013 and 2014 to coincide with the periods of grain filling where is the period where the ARI’s model is applied. In this research, the method of using a white reference panel described by McCoy KR and was adopted Milton EJ [9, 42]. The GER1500 was used to acquire measurements on the target and on the control panel. By applying the ratio of the reflected radiance from the target to the reflected radiance from the panel and by taking into account the control panel correction, the reflectance of the target was obtained. Field spectroradiometric measurements were collected in-situ using the GER-1500 field spectroradiometer with reflectance spectrum ranging from 350 nm to 1050 nm. The instrument’s field of view was 4o so the user was taking measurements of 1,2 m above the target (8 cm).

Click to enlarge

Using the GER1500 spectro-radiometer in situ, the surface reflectance values equivalent to the Landsat 8 OLI or/and ETM+ RED and NIR bands were determined. To filter the data through the Relative Spectral Response (RSR) values of Landsat 8 OLI/ETM+, the GER1500 reflectance values were interpolated to obtain the reflectance values at the incremental wavelength of the RSR. This was done since the GER1500 reflectance values were given at a different incremental wavelength scale. Then, the GER 1500 experimental data was filtered through the RSR functions and averaged within the limits of the RED and NIR 8 OLI/ETM+ bands, to yield the in- band reflectance values. Vegetation indices were determined for each date of acquisition and averaged to have the average value of the grain filling stage At each time, about 25 in situ measurements well spread in each plot, were collected in order to have a representative sample from the plot. Then the 25 measurements were averaged into a single representative measurement which was then processed to retrieve NDVI and WDVI. Simultaneously to spectroradiometric measurements, LAI measurements were taken using the SUN-SCAN canopy analyzer (Delta-T Devices Ltd., UK). The same procedure for averaging the LAI data has been followed and averaged LAI values have been created for each of the plots for the grain filling stage for each year. b) The field data was collected on schedule based on the satellite overpass in order to have comparable to the satellite imagery data. The satellite overpass was at the specific dates for 2013 (5, 13 and 21 of April for ETM+, 8 OLI and ETM+ respectively) and for 2014 (1 and 16 of April for 8 OLI and 8 OLI respectively). Thus averaged spectroradiometric data and averaged Landsat’s remotely sensed data regarding NDVI and WDVI of the grain filling stage of the Durum wheat have been prepared for each of the 25 plots in the area of interest. The WDVI index has the advantage to reduce to a great extent the influence of soil background on the surface reflectance values. Although simple, WDVI is as efficient as most of the slope based VI. The effect of weighting the red band with the slope of the soil line is the maximization of the vegetation signal in the near-infrared band and the minimization of the effect of soil brightness. After creating a set of data from bare soil of the area during the year, the slope of the line was set to 1.27. The slope line is created from the NIR and R spectrum using spectroradiometric data (200 sample random measurements) from the area’s soil (Figure 4).

Click to enlarge

ERDAS imagine (version 11) was used for the pre-processing and post-processing of satellite images. All satellite images were geo-referenced using a second order polynomial with 12 ground control points. Then all images have undergone atmospheric correction using the Darkest Pixel algorithm. The algorithm of Darkest pixel has been applied to each satellite images as for atmospheric effects removal. After the Digital Numbers (DN) of the selected dark target (Asprokremmos dam in the area of interest) near the area of interest are converted to units of radiance using calibration offset and gain parameters, the target reflectance at ground level is determined using the following simplified equation:

• Where
• $\rho g$ is the target reflectance at the ground
• Lds is the dark object radiance at the sensor
• Lts is the target radiance at the sensor,
• Eo.d is the solar irradiance at the top of the atmosphere corrected for earth-sun distance variation
• $\theta$ is the solar zenith angle

It is important to mention that dark object radiance at the sensor corresponds to the atmospheric path radiance which is subtracted or added from each spectral band from the target radiance at the sensor.

Also, LAI of the specific plots was estimated to be compared to the real values LAI which were measured from field measurements. For estimating LAI the following equation was applied [4; 43]:

$$LAI = -\frac{1}{a} \ln \left( 1 - \frac{WDVI}{\rho \infty (\lambda nir)} \right)$$

• Where
• LAI = Leaf Area Index,
• WDVI = Weighted Difference Vegetation Index,
• $\alpha$ = complex combination of extinction and scattering coefficients,
• $\rho \infty$ = asymptotically limiting value of the WDVI at very high LAI values.

Standard Values for $\alpha$ and $\rho \infty$ were taken from literature specifically for cereals, and used for the model [44;45].

The specific equation was applied to satellite images and LAI maps were created. Various studies have shown that empirical equations of LAI and Vegetation Indices can provide good estimates of LAI from the satellite images [4, 37, 43, 44, 45, 46].

c) Then the Crop Yield Prediction model has been applied by using atmospheric corrected and non-corrected available satellite imagery of April 3013 and 2014 mentioned earlier. In the model application LAI and NDVI are estimated using the satellite images while the other two factors, namely soil humidity and crop temperature were directly measured and averaged in each plot. Leaf Area Index (LAI), crop temperature and soil humidity at 40cm in situ measurements were also taken simultaneously to spectroradiometric measurements. Using the SUN-SCAN canopy analyzer (Delta-T Devices Ltd., UK), users can acquire the LAI value directly by setting it up perpendicular to the crop rows [47]. As well, a calibration measurement above the canopy is needed for each time the user takes a measurement in the crops. The SunScan was calibrated under a standard light source against an accurate Photosynthetic Active Radiation (PAR) quantum sensor. The spectral and cosine responses of the sensors approximate to the ideal response, but fall off at the extremes of the range. Under most normal daylight conditions, errors due to the deviation are small, but it is possible, in conditions such as under artificial light, to find larger errors in the absolute values measured. Since the Sunshine Sensor and Probe are closely matched, this has minimal effect on the canopy calculations which are based on ratios of incident and transmitted light [48]. Crop temperature (CT) was derived from a digital handheld thermometer while soil humidity (SH) was derived from soil moisture WET-2 from Delta Devices. Then, values of real Crop Yield along with estimated Crop Yield from corrected and not corrected satellite images have been created for every single plot for each year (2013 and 2014). These values, again, are the average values for the specific period of grain filling. The final product of the model application is a Crop Yield map for each satellite image, where values of Durum wheat yield are retrieved and recorded in the tables next to the real yield values, given from the Agricultural Research Institute of Cyprus.

d) Finally for each factor affecting Crop Yield but also for Crop Yield, a statistical method for comparing the paired values have been employed, in order to test if there is significant differences among the different pairs. Student’s t-test is used to find out basically if there is statistically significant difference between real and corrected values, real and uncorrected values and lastly between corrected and uncorrected values. The Student’s t-test for paired samples, as a statistical significance test, was applied to compare and assess the results meaning. To assess the value of t, the standard deviation SD of each pair of values (real and predicted) must be known. Dj refers to the difference of each pair, Da is the average difference of each pair and n is the number of the pairs.

2 D S Variance ( ) is calculated from equation (1):

(𝜮𝑫𝒋)𝟐

𝟐−

𝟐= 𝜮𝑫𝒋

𝒏 / 𝒏−𝟏) (1) 𝒔𝑫

2 Da S Following, the average Variance ( ) of all pairs must be calculated from the following equation (2)

𝟐 = 𝒔𝑫

𝟐/𝒏 (2) Finally, t value is calculated solving equation 3:

𝒔𝑫𝒂 𝑫𝒂

𝒕= 𝑺𝑫𝒂 (3)

Vegetation Indices and LAI comparison

Vegetation indices were calculated using the spectrum range acquired from the spectroradiometric measurements. Figures 5, 6 and 7 show the values of radiometric, uncorrected and corrected from atmospheric effects NDVI, WDVI and LAI respectively.

In Figure 5, the reference value is the NDVI radiometric which are basically the real values of NDVI.

Standard Deviations taken from Table 4, for 2013 and 2014, are applied in the figure in the form of discontinue line (perpendicular blue line) above and below the values. This is done to configure the movement of the corresponding Uncorrected (red line) and Corrected values (green line). It is very obvious that NDVI uncorrected values or line do not follow the schematic sequence of the reference values (radiometric values) but has a rather significant different movement from point to point. Most of the cases are even out of the limit of the Standard deviation of real NDVI meaning that there is a basic difference. On the other hand, NDVI corrected values or line has almost identical way as the reference line. All the cases are included in the limits of real NDVI Standard Deviation identifying that there are no significant differences between real and uncorrected values for NDVI. This is the case for years, 2013 and 2014. On the other hand values of uncorrected NDVI seem to have significant differences from the real (radiometric) but also the corrected values of the NDVI. Most of the points are out of the Standard Deviation limits implying that if there is a relation this could not be sound.

Click to enlarge

The same results come to light for the other two factors, WDVI and LAI. Uncorrected values of WDVI and LAI are differentiated from the reference (real) values or lines of the corresponding factors significantly and are out of the Standard Deviation’s limits for each factor, most of the times for the two years.

Click to enlarge

Corrected values are following the reference lines and are included in the Standard Deviation’s limits of each factor, except few single cases of WDVI 2014 and LAI 2013. Thus the general idea one receives from the specific figures and analysis is that uncorrected from atmospheric effects values have a significant difference from the real values while the corrected from atmospheric effects values are very close and within the limits of the real values’ Standard Deviation.

Click to enlarge

Figures 5-7 have provided a clue of the assumption raised by the authors that atmospheric effects have negative effect on the models using remotely sensed data. It has been shown that values from satellite images that have not undergone any atmospheric correction will have significant difference from the real values. Now intended purpose is to verify that statistically, and to illustrate that corrected from atmospheric effects values have no significant differences from the real values.

Student’s t-test was applied based on the results (paired values) of the NDVI, WDVI and LAI for the two years. The results of the test are shown in Table 5.

SPSS statistical software was used to proceed with the statistical analysis and obtain the values of T-test to be compared with those of existing Statistical Tables for T- test shown at the third column of Table 5. The analysis has illustrated that the values of T observed between real values and corrected values of NDVI, WDVI and LAI were smaller than the Tstatistical, which implies that for (n-1) degrees of freedom and at a confidence level of 95%, these values do not have a significant statistical difference between them. Thus for all the cases where satellite images are atmospherically rectified there is not a statistically significant difference between the estimated (corrected) and measured (real) values.

In the case of comparing real NDVI, WDVI and LAI values to uncorrected values, the opposite phenomenon takes place, meaning that there is, in all cases, a statistically significant difference of the paired values, which means that uncorrected images cannot be used to After applying the Drum wheat prediction model yield values for 2013 and 2014 were retrieved. These values were compared graphically and statistically to real values in order to test the null hypothesis that if satellite images are not compensated for atmospheric correction then they cannot be used in crop yield modelling.

retrieve the corresponding NDVI, WDVI and LAI values since they will be misleading.

Crop Yield comparison

Click to enlarge

Figure 6, illustrates the real values of Crop Yield given from the Agricultural Research Institute of Cyprus for the specific plots which are considered to be the reference values. Standard Deviation of the specific values taken from Table 4, for 2013 and 2014 and applied in the figure in the form of discontinue line (perpendicular blue line) above and below the values.

This is done to configure the movement of the corresponding Uncorrected (green line) and Corrected values (red line). It is very obvious that Crop Yield uncorrected values or line do not follow the schematic sequence of the reference values (real values) but has a rather significant different movement from point to point. Most of the cases are even out of the limit of the Standard deviation of real values meaning that there is a basic difference. On the other hand Crop Yield corrected values or line has almost identical way as the reference line. All the cases are included in the limits of real values’ Standard Deviation identifying that there is no significant differences between real and uncorrected values. This is the case for years, 2013 and 2014.

Finally, Student’s t-test was applied based on the results (paired values) of the NDVI, WDVI and LAI for the two years. The results of the test are shown in Table 6.

The analysis has illustrated that the values of T observed between real values and corrected values of Crop Yield were smaller than the Tstatistical, which implies that for (n- 1) degrees of freedom and at a confidence level of 95%, these values do not have a significant statistical difference between them. Thus for all the cases where satellite images are atmospherically rectified there is not a statistically significant difference between the estimated (corrected) and measured (real) values.

In the case of comparing real Crop Yield values to uncorrected values, the opposite phenomenon takes place, meaning that there is, in all cases, a statistically significant difference of the paired values, which means that uncorrected images cannot be used to retrieve the corresponding NDVI, WDVI and LAI values since they will be misleading.