Assessment of Maize Growth and Development with High- and Medium-Resolution Remote Sensing Products

Abstract

The availability of a great amount of remote sensing data for precision agriculture purposes has set the question of which resolution and indices, derived from satellites or unmanned aerial vehicles (UAVs), offer the most accurate results to characterize vegetation. This study focused on assessing, comparing, and discussing the performances and limitations of satellite and UAV-based imagery in terms of canopy development, i.e., the leaf area index (LAI), and yield, i.e., the dry aboveground biomass (DAGB), for maize. Three commercial maize fields were studied over four seasons to obtain the LAI and DAGB. The normalized difference vegetation index (NDVI) and visible atmospherically resistant index (VARI) from satellite platforms (Landsat 5TM, 7 ETM+, 8OLI, and Sentinel 2A MSI) and the VARI and green canopy cover (GCC) from UAV imagery were compared. The remote sensing predictors in addition to the growing degree days (GDD) were assessed to estimate the LAI and DAGB using multilinear regression models (MRMs). For LAI estimation, better adjustments were obtained when predictors from the UAV platform were considered. The DAGB estimation revealed similar adjustments for both platforms, although the Landsat imagery offered slightly better adjustments. The results obtained in this study demonstrate the advantage of remote sensing platforms as a useful tool to estimate essential agronomic features.

1. Introduction

Maize (Zea mays L.) is the most produced cereal in the world, with more than one thousand million tons harvested in 2019 [1]. The maize harvested area was close to 200 million hectares, representing more than 14% of the worldwide cropped land. It is the most important crop for grain production, although wheat and rice are the most important for direct human consumption [2]. Research has forecasted that the world population will grow to reach 9.1 billion in 2050. Therefore, maize should be considered one of the main pillars for global food security.

It is not possible to speak about global population increase without considering future global warming scenarios where the water availability will decrease in amount and quality. Maize water requirements range from 500 to more than 800 mm year⁻¹ [2], which makes this cereal one of the most water demanding. Most maize cultivated lands require supplemental irrigation for optimal growth. Therefore, maximizing the crop production and reducing the water applied for the same cultivated land (crop water productivity) is one of the main challenges for global agriculture. In this context, precision agriculture plays a key role by improving sustainable decision making, defined as a way to apply the right treatment in the right place and at the right time [3].

Precision agriculture demands crop monitoring and management as well as crop biomass and yield estimation, considering the leaf area index (LAI, m² leaf m⁻² soil) as one of the most representative indicators [4]. With the objective of improving yield predictions, accurate LAI characterization constitutes one of the first steps for crop modeling and, therefore, helps farming decision making. LAI characterization along the growth cycle provides significant information regarding the canopy architecture, efficiency of light interception, and yield—primarily for those crops such as maize with a large photosynthetic system. Field measurements are commonly destructive and tedious, and sampling is not repeatable across time. Intensive work has been conducted to obtain non-destructive approaches based on the empirical relationships between LAI and spectral vegetation indices (VIs) [5] or crop structure parameters, such as green canopy cover (GCC) (Ballesteros et al. 2014).

Satellite remote sensing is widely used for crop monitoring through the relationship between the normalized difference vegetation index (NDVI) and crop vigor [6], where infrared bands are essential for NDVI calculations. High-spatial-resolution observations are not obtained with conventional satellite missions, such as Landsat or Sentinel. A new generation of satellites, e.g., RapidEye or Formosat-2, offer higher revisit cycles and resolutions, although they do not supply open data. Sentinel-2 appeared as an opportunity to improve these concerns with the effectiveness of decametric resolution satellite imagery to describe and monitor crop growth and development [6]. This open data satellite imagery is not suitable for field monitoring purposes in agricultural systems as it causes mixed soil, crop, and weed pixels, leading to biased VIs [7]. Using different platforms and sensors demands the analysis of interoperability between the different products. The authors of [8] compared the atmospherically corrected NDVI from Landsat 8 and that from Sentinel-2. They found that Sentinel-2 offered slightly lower values than Landsat 8 did, with a very high correlation coefficient. However, the RMSE was 0.06, which is lower than the generated model errors. We can therefore be confident that this aspect did not highly impact the final results. The NDVI was developed to try to reduce these negative impacts. The interoperability between Landsat 7 and 8 [9] showed no deviation for high NDVI values and slight deviation when sparse vegetation was found (early crop development stages). Additionally, [10] showed an RMSE of 0.08 for the NDVI TOA between Landsat 8 and Landsat 7. Again, we found a low impact of this deviation compared with the accuracy of the final model. When analyzing the interoperability between Landsat 7 and Landsat 5, [11] indicated a high degree of similarity, which implies that monitoring activities initiated using Landsat 5 data can be continued with a minimal amount of caution using Landsat 7 data. In the case of UAV-based information, vegetation segmentation was not performed. The spectral response of soil and vegetation was averaged for each sampling area. We followed this strategy to ensure interoperability between satellite-based and UAV-based data.

The studies of Haboudane et al. (2004) and Viña et al. (2011) are only two examples of the great efforts made to describe the LAI using spectral data since remote sensing technologies emerged in the early 1980s. VIs have demonstrated a high sensitivity to the LAI. The NDVI has been one of the most traditionally used indices. It is very well known, as the soil background influences LAI characterization when using NDVI [12,13]. Other authors utilized other indices, such as the soil-adjusted vegetation index (SAVI), to overcome this problem [14].

The visible atmospherically resistant index (VARI) has been widely described in the scientific literature for crop monitoring through different VIs, such as LAI estimation [15]. In the case of satellite-based remote sensing, where low spatial resolution appears, VARI is not widely utilized because it is based on the relationship between visible bands, and the response of vegetation in the visible spectrum is much lower than that in the infrared region of the spectrum.

However, in the case where a higher resolution is obtained, such as in the case of UAV-based remote sensing, the visible spectrum supplies useful information that is not contaminated by elements other than vegetation [16]. However, Haboudane et al. agreed that it is difficult to recommend a specific index that describes only the desired variable and is insensitive to the other environmental or biological parameters [12].

Unmanned aerial vehicles (UAVs), are becoming a pillar for precision agriculture purposes in the context of Agriculture 4.0. They are easy to use, acquiring a platform is becoming more and more affordable, the obtained products are independent of cloudy conditions, and they are able to capture images of a higher spatial and temporal resolution [17] compared with the traditionally used Landsat or modern Sentinel satellites. However, there are still some important limitations, e.g., in large-scale applications. Large-scale remote sensing data allow for monitoring the heterogeneity of cropland, whilst high temporal resolutions allow for accurate description of crop growth, especially when farming management decisions are required [4].

Most of the studies based on information derived from UAVs utilize multispectral information, which is an attempt to adapt the indices and methodologies used for satellite-based remote sensing to the images captured by UAVs. The drastic change in the spatial resolution must be accompanied by a change in the approach to obtain useful information for agriculture. The choice of red, green, and blue (RGB) bands when using UAV platforms is one of the main changes, with more affordable sensors and easier photogrammetry processes. The use of RGB information can be exploited by generating new RGB-based indices [15,16] or even exploiting the generated point clouds and, therefore, the 3D information [16,18]. Many scientific works achieved significant results by using RGB cameras mounted on UAV platforms [19,20,21,22].

Currently, real-time crop growth and biomass prediction, where LAI plays a key role [23], is one of the main concerns not only for research purposes but also for farming decision management. Many research works have been developed with this aim [5,6,21,22,24]. Most of them focused on studying specific resolutions, high or medium (satellite platforms) and high or very high resolutions (UAV platforms). However, to the authors’ knowledge, comparisons of both platforms have not been sufficiently developed for most agricultural frameworks, as most of them were focused on vineyards, orchard agricultural systems, and forests [7,25].

The need to answer the question of whether a higher resolution is needed in addition to deciding which indices—derived from satellite or UAVs—offer more accurate results to characterize maize vegetation inspired this present research. The objective of this research work was to assess, compare, and discuss the performance and limitations of satellite-based and UAV imagery in terms of the canopy development (LAI) and dry aboveground biomass (DAGB) for maize. To achieve this objective, a reliable description of the crops is required.

2. Materials and Methods

2.1. Study Area

This study was carried out in Tarazona de La Mancha (Albacete) in the southeast of Spain (Figure 1) over four seasons in 2011, 2012, 2015, and 2016. This area is classified as semi-arid, as the aridity index reaches 0.26 [22]. The meteorological variables were recorded from an agro-meteorological station located close to the field sites (X = 593,176 m; Y = 4,346,196 m (ETRS89, UTM30N EPSG: 25830), 722 m asl).

The daily minimum temperature (T_min) values for the maize growing season ranged from 7.7 to 9.8 °C in 2011 and 2016, respectively; the daily maximum temperatures (T_MAX) were higher than 33 °C (Table 1). According to registered agro-climatic data, the 2016 season was dry compared to the rest of the seasons, with precipitation (P) and reference evapotranspiration (ETo) values higher than the main registered ones for the four seasons (P = 230 mm and ETo = 942 mm) (Table 1).

Field experiments were performed in three different fields of irrigated FAO-700 cycles of maize grain crop: maize was grown in field A (18.02 ha) for the 2011 and 2012 seasons, in field B (42 ha) for the 2015 season, and in field C (37 ha) for the 2016 season (Figure 1). Field A was cultivated under conventional farming practices, whilst fields B and C were cultivated following conservation agriculture practices, i.e., promoting minimum mechanical and soil disturbances and permanent soil organic coverages. The landscape is characterized by a vast plain. The soil is classified as Inceptisol Xerept in the study area [26,27]. Soil analyses were performed before seeding for every season. The soil was characterized as sandy clay loam for field A and as sandy loam for fields B and C. Local soils were characterized by alkaline pH values, i.e., higher than 8.0. Organic matter contents ranged from 0.91 to 1.58% for field A and from 3.8 to 4.2% for fields B and C, respectively. The higher values of organic matter content for fields B and C compared with field A show how no tillage systems increase organic matter content. The maize was sown on 29 April 2011, 26 April 2012, 3 April 2015, and 29 March 2016 with a seed density of 120,000 seeds ha⁻¹ for the four seasons. The field experiments were performed in the mentioned three commercial fields belonging to different landowners under different cultivation techniques.

The growing degree days (GDD) were calculated using the double sine wave method [28,29]. The minimum and maximum threshold temperatures were 6 and 26 °C for the FAO-I stage, 7 and 31 °C for the FAO-II stage, 12 and 38 °C for the FAO-III stage, and 6 and 40 °C for the FAO-IV stage, respectively [30,31]. Harvest occurred on 28 and 30 November in 2011 and 2012 and 18 and 28 September in 2015 and 2016.

Although the harvest dates were quite different for 2011 and 2012 compared with 2015 and 2016, the four seasons reached full ripeness when more than 2000 GDD were accumulated. Field A was irrigated using a permanent solid system (4918 m³ ha⁻¹ in 2011 and 5210 m³ ha⁻¹ in 2012). Fields B and C were irrigated using a central pivot system (8816 m³ ha⁻¹ and 8067 m³ ha⁻¹). Differences in the total amounts of irrigation were related to different irrigation strategies and scheduling as well as different farm managing strategies. These differences in farming strategies will serve to evaluate the generalizability of the obtained relationships.

2.2. Field Data Collection

For the four analyzed seasons, the maize was monitored weekly according to the Bundesanstalt Bundessortenamt and Chemical Industry scale (BBCH) [32]. The yield values were 15,300, 13,900, 16,674, and 17,343 kg ha⁻¹ for the 2011, 2012, 2015, and 2016 seasons, respectively, according to the last sampling event when full ripening was reached.

The dates for field data acquisition were selected according to the main canopy changes and considering the main phenological stages as defined by [33]. The dates, together with the main phenological stages, are presented in Table 3. The sampling plots were randomly selected for every date. Coordinates are shown in Supplementary Table S1. Three 1.4-square-meter sampling frames delimited each sampling plot for the 2011, 2012, and 2016 seasons and four for the 2015 season. The total maize plants per frame were counted to determine the plant density. Six plants per sampling plot were collected. The LAI was calculated, determining the leaf area (LA) using an automated infrared imaging system, LI-COR-3100C (LI-COR Inc., Lincoln, NE, USA), and the plant density (Figure 2, step 3). The sampled maize plants were carried to the laboratory where the dry aboveground biomass (DAGB) was determined for every sampling event. The coordinates of the sample points were obtained with a GNSS-RTK with an accuracy of 2 cm in planimetry.

2.3. UAV Imagery

Very-high-resolution flights and field sampling events were performed at the same time. Red, green, and blue band (RGB) aerial images were taken using a Microdrone md4-200 (Microdrones, Inc., Kreuztal Germany) provided with an RGB PENTAX A40 digital camera (Pentax, Golden, CO, USA) in 2011 and 2012 and with a Microdrone md4-1000 (Microdrones, Inc., Kreuztal, Germany) provided with a SONY α ILCE-5100L (SONY Corporation, Tokyo, Japan) in 2015 and 2016. To obtain a ground resolution of 1.03 cm pixel⁻¹, flights were performed at a 44-meter height for the md4-200 and 48 m for the md4-1000 microdrone.

Aerial images were taken automatically following a flight plan for each UAV platform following the Microdrones Photogrammetric Flight Planning software (MFLIP) output [34] (Figure 2, step 1.2.). This software requires, as input parameters, a digital elevation model (DEM) of the area of flight and the internal orientation parameters of the digital camera. This software also considers GPS errors, which results in an improvement of the flight planning performance. Proper overlapping imagery was established in 60% and sidelapping was established in 25%. Blurred images were automatically detected following the proposed methodology of [35]. The geomatic products (orthoimage, digital surface model (DSM), and point cloud) were obtained using PhotoScan Professional version 1.4.1 software (Agisoft LLC, St. Petersburg, Russia). Flights were always performed at solar midday and sun glints and hotspots were detected following the advice of [36].

2.4. Satellite Data Acquisition and Preprocessing

Satellite images were downloaded considering the closest dates to sampling events and UAV flights for every season, as climatic conditions, such as clouds, did not allow for obtaining low and medium imagery (satellite) from when sampling events and UAV flights were performed for the same dates. The analyzed satellite images were extracted from operational Landsat 5 TM (2011 season), Landsat 7 ETM+ (2011 and 2012 seasons), Landsat 8 OLI (2015 and 2016 seasons), and Sentinel 2A MSI (2016 season) (Figure 2, step 1.1.). Reflectance images with atmospheric correction were used. The Landsat images were download free of charge from the United States Geological Survey (USGS) [37]. Landsat 8 OLI, Landsat 5 TM, and Landsat 7 ETM+ collection 2 surface reflectances are Level-2 data products. Level-2 data products present radiometric calibration and atmospheric correction algorithms to Level-1 Landsat data products. In the case of Landsat 8 OLI, the data were processed by the suppliers using the Land Surface Reflectance Code (LaSRC). Landsat 5 TM and Landsat 7 ETM+ data were processed by the suppliers using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [37]. Sentinel images were from the Copernicus Open Access Hub website [38] and corresponded with the Level-2A product (orthorectified Bottom-Of-Atmosphere reflectances) The Level-2A product presents an atmospheric correction applied to Top-Of-Atmosphere (TOA) Level-1C orthoimage products [39]. The Level-2A Sentinel 2 products were processed by the suppliers using Sen2Corprocessor. The final products present geometric, radiometric, and atmospheric corrections. Specifically, 43 cloud-free images were used in this work (Table 2).

2.5. Calculation of the Vegetation Indices and Parameters

In the case of satellite-based remote sensing, the relationships between the LAI and DAGB field data and the VIs, NDVI_Landsat, NDVI_Sentinel, VARI_Landsat, and VARI_Sentinel were computed. The NDVI is the most widespread VI, while VARI was calculated to compare the results of using the RGB bands obtained with the drones with the visible bands for each utilized mission. The NDVI was computed using Formula (1) [40], where NIR and R are the near-infrared and red bands of the image, respectively (Figure 2, step 2.1.).

(1)$NDVI=\frac{\left(NIR-R\right)}{\left(NIR+R\right)}$

VARI was determined following Gitelson et al. (2002) (2), where r, g, and b are the normalized red, green, and blue bands of the image, respectively (Figure 2, step 2.2.).

(2)$VARI=\frac{\left(g-r\right)}{\left(g+r-b\right)}.$

The NDVI and VARI values were linearly interpolated between adjacent dates (Table 2) to obtain the daily values for every field for the same dates on which sampling events and UAV flights were performed (Table 3) analyzed in the study.

In the case of UAV-based remote sensing data, VARI_UAV and GCC_UAV were calculated to monitor the LAI and DAGB along the growth cycle (Figure 2, step 2.2.). The radiometric quality of conventional RGB cameras is much lower than the radiometric quality provided by satellite missions. Thus, to minimize this effect, before calculating the proposed VI, digital numbers (DNs) were calculated from the average of the digital number of the R, G, and B bands of every sampling plot. Nevertheless, non-normalized RGB coordinates are directly proportional to the total light reflected from a surface and are highly sensitive to the intensity of the illuminating source and its angle with the target surface. To minimize this problem, the RGB coordinates (the DNs for every band) were used to derive the red, green, and blue chromatic coordinates (ρ, γ, and β) [16,41], ensuring normalized values for each band. This procedure also minimizes the negative effect of using different sensors for the different years. The GCC was determined using the Leaf Area Index Calculation software [42], which segments the vegetation in an orthoimage using a machine learning algorithm.

2.6. Statistical Analysis

Statistical analysis was performed using MATLAB™ (MathWorks, Natick, MA, USA) (Figure 2, step 4). Two different groups of predictor variables were used to assess the LAI and DAGB: (i) a predictor derived from the agro-climatic variable, i.e., GDD (Figure 2, step 2.1. and step 2.2.); and (ii) the previously mentioned VIs, i.e., NDVI and VARI for the satellite imagery and GCC and VARI for the UAV imagery (Figure 2, step 2.1. and step 2.2.).

A multilinear regression model (MRM) was explored to estimate both target variables, i.e., LAI and DAGB, for the 2011, 2012, and 2015 seasons from Landsat missions (Landsat 5 TM, Landsat 7 ETM+, Landsat 8 OLI) (NDVI_Landsat and VARI_Landsat) and the calculated GDD. The MRM was used to establish the relationships between the GCC_UAV and VARI_UAV obtained from the UAV imagery and the calculated GDD from the field data. As performed for the satellite approach, LAI and DAGB were trained using GCC_UAV, VARI_UAV, and GDD as independent variables for the 2011, 2012, and 2015 seasons. The testing of the selected models was performed using independent data from the 2016 season, which was not used in the previous processes.

The Sentinel platform was launched in 2015. Stepwise linear regression (SLR) was explored to predict the LAI and DAGB for the 2016 season to evaluate the performance of this mission. SLR may be considered when a small dataset is used, without losing a significant portion of the explanatory power of the data [43]. This was the case for the evaluation of the Sentinel 2A MSI mission, for which only one irrigation season (2016) was available in this study.

The considered predictor variables were NDVI_Sentinel, VARI_Sentinel, and the computed GDD for the LAI and DAGB. SLR was performed using forward and backward regression to determine the final model. The selected criterion to add or remove terms was the p-value for the F test. If the p-value of the model was smaller than 0.05, then the term was added to the model, whereas if the p-value of the model was higher than 0.1, then the term was removed from the model.

The adjusted coefficient of determination (R²_adj) and the root mean squared error (RMSE) were used to assess the performance of the obtained models.

3. Results

3.1. Field Data Sampling

The predictor variables NDVI, VARI, GCC, and GDD and the field sampling data LAI and DAGB for every sampling event are shown in Table 3. VARI_Landsat, VARI_Sentinel, and VARI_UAV for every sampling event in every season are shown in Supplementary Figure S1. LAI values ranged from 0.31 to 5.36 m² leaf m⁻². The maximum LAI values were sampled at the development of the fruit and ripening phenological stages, when maximum predictor variable values for NDVI_Landsat, NDVI_Sentinel, VARI_Landsat, VARI_Sentinel, GCC_UAV, and VARI_UAV were also reached (at around 1100 GDD). For the highest LAI values, the coefficients of variation (CVs) were the lowest.

The maximum DAGB reached 3473, 3479, 2910, and 2924 g m⁻² soil in 2011, 2012, 2015, and 2016, respectively. The maximum DAGB values were registered at ripening when more than 3400 GDD was reached for the 2011 and 2012 seasons and at around 2900 GDD for the 2015 and 2016 seasons. The 2012 season showed higher heterogeneity regarding the DAGB compared with the others: the CV reached 36% in 2012, while it went from 8 to 11% for 2011, 2015, and 2016. The maximum values of all predictors from the satellite and UAV platforms appeared at the same time that the LAI did: the third FAO stage. The NDVI_Landsat and NDVI_Sentinel maximum values ranged from 0.69 to 0.89. The VARI_Landsat and VARI_Sentinel maximum values reached 0.25. The GCC reached 95.36%, and the VARI_UAV reached 0.09.

3.2. Relationships between the Predictors, LAI, and DAGB

Table 4 shows the studied models’ estimation of the LAI from NDVI_Landsat and GDD, VARI_Landsat and GDD, GCC_UAV and GDD, and VARI_UAV and GDD for the 2011, 2012, and 2015 seasons. Better adjustments were obtained when the predictor variables from the UAV platform were considered. The GCC_UAV and GDD predictors offered values of 0.931 and 0.468 m² leaf m⁻² soil for the R²_adj and RMSE, respectively, whilst predictors from Landsat were lower than 0.730 for the R²_adj, and the RMSE values were higher than 1.00 m² leaf m⁻² soil.

DAGB estimation revealed similar adjustments for both platforms, although the Landsat imagery offered slightly better adjustments: The VARI_Landsat and GDD predictors showed values of 0.779 and 578 g m⁻² soil for the R²_adj and RMSE, respectively, whilst the NDVI_Landsat and GDD showed similar adjustments (0.770 and 589 g m⁻² soil for R²_adj and RMSE, respectively). As adjustments obtained for the Landsat platform, the DAGB model obtained by GCC_UAV and GDD offered values of 0.747 and 615 g m⁻² soil for the R²adj and RMSE, respectively, and the model described by the VARI_UAV and GDD predictors offered similar statistical adjustments compared with the GCC_UAV and GDD.

Figure 3 shows, graphically, the results obtained when the predicted and observed LAI values are represented for the 2016 season using GCC_UAV and GDD as predictor variables and the selected model. We observed that the LAI_predicted values were overestimated when they were lower than 3.48 m² m⁻² soil.

The predicted versus observed DAGB values using VARI_Landsat and GDD as predictor variables can be observed graphically in Figure 3 for the 2016 season. The statistical indicators revealed significant adjustments: R²adj = 0.875. The maximum deviation was observed at the end of the growing cycle when the crop was completely developed according to the sampling data (Table 4). An overestimation for the DAGB predicted values can be observed along the whole growing cycle using VARI_Landsat and NDVI_Landsat (Figure 4).

3.3. Relationships between LAI and DAGB and the Predictor Variables Obtained from Sentinel

Table 5 shows the obtained models for LAI and DAGB from the Sentinel imagery in 2016. LAI prediction by NDVI_Sentinel and GDD showed values of 0.871 and 0.540 m² m⁻² for the R²_adj and RMSE, respectively. The DAGB model revealed that the only significant predictor was GDD. SLR did not consider any of the VIs proposed, and only GDD was selected as a predictor. Therefore, considering NDVI_Sentinel or VARI_Sentinel as predictor variables offered the same model and no statistical differences.

4. Discussion

The results obtained in this study demonstrate the advantage of remote sensing platforms as a tool to estimate the essential agronomic features of maize fields under different management practices. The significant correlations found in this study between the VIs derived from remote sensing platforms and the agronomic features, such as LAI and DAGB, offer promising results in the context of climate change and its implications for food security.

High canopy uniformity was observed when the LAI was at the maximum, as the CV values were low when most plants were at their maximum development. Therefore, the heterogeneity regarding canopy development in the maize crop rates disappeared when the fruit development phenological stage was reached (full canopy reached). We observed higher heterogeneity for the DAGB compared with the LAI. The 2012 season offered the highest CV (36%), whilst this heterogeneity did not appear for the LAI variable (CV was 10%) for the same season.

These results agree with those obtained by [44], who observed high heterogeneity rates—higher than 20% CV when DAGB was measured. Higher rates of CV for the DAGB were noted at the end of the growing cycle when the maximum values were reached in the 2012 season, when the amount of water for irrigation was much lower than that applied in 2015 and 2016, although the ETo values were similar.

The VARI reached quite different values depending on the platform used when the maximum canopy development occurred. We observed that the VARI_Landsat and VARI_Sentinel offered significantly higher values than VARI_UAV did. It is important to highlight that the different ranges of VARI values extracted from the different platforms did not affect the goodness of the adjustments, as they were considered independent predictors for LAI and DAGB for every season, and they were not combined. VARI_Landsat appeared to be the best indicator for DAGB estimation. Additionally, VARI_UAV offered significant results. We observed the high dependency of the full canopy for these indices, as the CV was lower when the first growth and development stages were completed. The NDVI_Landsat and NDVI_Sentinel revealed almost no variations when the full canopy was reached, indicating a saturation of this index when the canopy closes [43,44].

The obtained results showed better adjustments between the UAV predictor variables than those obtained from the satellite to represent canopy development. The NDVI from satellite and crop canopy radiometric information can be altered by other sources, such as inter-row paths [6]. The generated information from satellite platforms could be insufficient to properly evaluate the canopy behavior. Therefore, the predicted LAI from the satellite VIs led to a bias at the early crop growth stages, when the influence of the background could bias the LAI predictions [25]. In this case, imagery with a very high resolution is required to properly assess the canopy development.

Discriminating soil from crops is feasible using raw RGB images, but quantifying the percentage of canopy cover is necessary to classify these classes. As shown in the results, GCC_UAV was highly sensitive in monitoring the canopy development. This predictor variable could potentially be used not only to determine the LAI but also to represent canopy development for irrigation scheduling and early decision farming management. As mentioned before, at the end of the growing cycle, when reaching the maximum values of LAI, the NDVI might not reflect the increase in the canopy density linearly [44]. In these cases, imagery with a very high resolution is required to properly assess the canopy development.

Despite the good results shown for the VARI_Landsat and GDD for DAGB prediction, this is not the only effective index for this estimation. NDVI_Landsat offered very good adjustments when it was combined with the GDD predictor. The differences between both indices were not highly significant. Therefore, we recommend the use of NDVI_Landsat and GDD to predict the DAGB, as has been traditionally used by remote sensing users.

The results obtained for training and validation of the DAGB model from the Sentinel data using SLR indicated that the DAGB was only determined by the GDD behavior. This suggests that the remote sensing information did not provide any information for DAGB modeling. The absence of significant relationships could be related to the need for more data to create a proper DAGB model, as was observed for the Landsat satellites.

Therefore, a multilinear regression model (MRM) was used to evaluate the potential of NDVI_Sentinel, VARISentinel, and GDD as predictors to calibrate a DAGB model (Table 6). As was described for Landsat imagery, both VIs, NDVI_Sentinel and VARI_Sentinel, combined with the GDD revealed excellent adjustments. Nevertheless, the obtained results showed the difficulty in producing an accurate model to predict the DAGB, suggesting the necessity of creating a DAGB-predicting model in real time. The model should be fed the captured remote sensing data and recorded GDD training to test the model for every season in real time.

5. Conclusions

Accurate prediction of the main agricultural features for maize is a pillar for agricultural management overall in those countries where this crop is essential for food and feeding. Remote sensing platforms offer a huge amount of data for farming decision support tools. Comparing the accuracy of the different available platforms, i.e., satellites and UAVs, and the different generated vegetation indices and parameters is mandatory for better canopy and yield characterization. The LAI can be successfully characterized by using the indices derived from UAV platforms in addition to phenological information (GDD).

Whilst the predictors from Landsat satellites and UAV platforms offered similar results, the obtained results revealed that more data are required to evaluate the potential of the Sentinel data for LAI and DAGB prediction. The difficulties in performing an accurate prediction of DAGB suggest the importance of DAGB in real time when feeding a model with captured and processed remote sensing data. Satellite platforms offer open data and cover vast areas, whilst UAVs are more adaptable to different weather and field conditions. The use of RGB bands instead of traditional multispectral ones derived from the use of satellite data is one of the main opportunities with more affordable sensors and easier photogrammetry processes. Applying this methodology, LAI and DAGB prediction models must be calibrated for many other crops in different irrigation areas.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/agronomy11050940/s1, Table S1: Sampling events coordinates (ETRS89, UTM30N EPSG 25830), Figure S1: Visual Atmospheric Resistance Index from unmanned aerial vehicle values for 2011, 2012, 2015, and 2016 seasons, Visual Atmospheric Resistance Index from Landsat values for 2011, 2012, and 2015 seasons, and Visual Atmospheric Resistance Index from Sentinel values for 2016 season.

Author Contributions

R.B. is the principal investigator and corresponding author. She led and supervised the overall research and field data collection. J.F.O. wrote the paper in collaboration with R.B. and M.A.M., and F.B. performed the field data collection and UAV flights. L.G.-G. processed the satellite imagery. M.A.M., R.B., and J.F.O. performed financial and contractual duties. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Castilla-La Mancha Regional Government SPLY/19/180501/0080 and FEDER funds.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Locations of the studied maize commercial fields and agro-climatic stations in Tarazona de La Mancha (Spain).

Figure 2. The workflow from the flight image acquisition to the statistical analysis. NDVILandsat: normalized difference vegetation index from Landsat imagery; VARILandsat: visible atmospherically resistant index from Landsat imagery; NDVISentinel: normalized difference vegetation index from Sentinel imagery; VARISentinel: visible atmospherically resistant index from Sentinel imagery; GDD: growing degree days; GCCUAV: green canopy cover from unmanned aerial vehicle imagery; VARIUAV: visible atmospherically resistant index from unmanned aerial vehicle imagery; LAI: leaf area index; DAGB: dry above ground biomass; MRM: multilinear regression model; SLR: stepwise linear regression; R2adj: adjusted coefficient of determination; RMSE: root mean squared error.

Figure 3. Results of the estimation of the leaf area index (LAI) for the testing irrigation season (2016).

Figure 4. Results of the estimation of the dry aboveground biomass (DAGB) for the testing irrigation season (2016) using DAGB = −174.24–1317.3.1NDVILansat + 2.3997GDD.

Table 1

The main registered agro-climatic data throughout the four seasons (2011, 2012, 2015, and 2016).

Season	Tmin	TMAX	P	ETo
2011	9.8	33.3	215.8	921.0
2012	9.2	33.8	311.6	899.9
2015	9.6	36.2	219.5	954.1
2016	7.7	33.6	175.3	995.5
Annual mean	9.1	34.2	230.6	942.6

T_min: mean daily minimum temperature (°C); T_MAX: mean daily maximum temperature (°C); P: total rainfall (mm); ETo: total reference evapotranspiration (mm).

Table 2

Dates of the images used in the analysis.

Field A (2011)		Field A (2012)		Field B (2015)		Field C (2016)
L5	9 April 2011	L7	10 April 2012	L8	6 May 2015	S2A	1 May 2016
L5	11 May 2011	L7	28 May 2012	L8	22 May 2015	S2A	21 May 2016
L5	19 June 2011	L7	13 June 2012	L8	7 June 2015	L8	9 June 2016
L5	28 June 2011	L7	8 July 2012	L8	30 June 2015	S2A	13 June 2016
L7	7 August 2011	L7	31 July 2012	L8	9 July 2015	S2A	20 June 2016
L7	23 August 2011	L7	25 August 2012	L8	16 July 2015	S2A	23 June 2016
L5	7 September 2011	L7	3 October 2012	L8	1 August 2015	L8	2 July 2016
				L8	26 August 2015	S2A	10 July 2016
				L8	11 September 2015	L8	11 July 2016
						L8	18 July 2016
						L8	27 July 2016
						S2A	2 August 2016
						S2A	9 August 2016
						L8	12 August 2016
						L8	19 August 2016
						S2A	22 August 2016
						S2A	29 August 2016
						L8	4 September 2016
						S2A	8 September 2016

L5: Landsat 5; L7: Landsat 7; L8: Landsat 8; S2A: Sentinel 2A.

Table 3

Statistics of the field sampling data, satellite vegetation indices, and unmanned aerial vehicle (UAV) indices and data.

Season	Sampling Date Event	Principal Growth Scale (BBCH Scale) a	GDD	Field Sampling Data				Satellite Vegetation Indices								UAV Indices and Data
				LAI		DAGB		NDVILandsat		VARILandsat		NDVISentinel		VARISentinel		GCC		VARIUAV
				$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV	$\overline{x}$	CV
2011	22 June	1: Leaf development	409	0.85	9.74	66.12	19.27	0.43	3.42	−0.13	−3.49					21.44	4.43	−0.23	−42.45
	20 August	7: Development of fruit	1244	5.14	6.07	996.30	20.97	0.69	0.00	0.07	0.00					88.25	1.98	0.09	19.85
	30 August	8: Ripening	1412	4.78	22.09	3473.17	8.41	0.61	0.00	−0.04	0.00					78.65	22.20	−0.03	−57.96
2012	21 June	1: Leaf development	428	1.54	0.00	136.07	0.00	0.40	0.00	−0.13	0.00					34.84	0.00	−0.11	0.00
	12 July	6: Flowering, anthesis	698	4.36	8.39	1333.65	27.05	0.65	0.00	−0.01	0.00					75.32	16.75	0.08	22.45
	14 August	7: Development of fruit	1198	4.61	5.69	2096.95	2.12	0.71	0.00	0.09	0.00					95.36	0.16	0.05	15.60
	28 August	8: Ripening	1455	3.84	18.90	2750.07	16.76	0.68	0.73	0.09	20.61					58.56	6.38	−0.03	−21.79
	10 September	8: Ripening	1640	2.88	29.09	3478.70	36.24	0.52	1.47	−0.04	−16.95					37.23	6.50	−0.09	−21.88
2015	29 May	1: Leaf development	403	0.44	14.46	32.33	20.40	0.34	66.69	−0.06	−66.83					12.60	7.19	−0.23	−2.81
	10 June	3: Stem elongation	576	2.14	9.41	191.33	10.12	0.60	11.24	0.02	90.36					47.98	8.65	−0.03	−47.38
	25 June	6: Flowering, anthesis	784	5.36	10.42	658.52	19.20	0.58	66.69	0.08	66.68					84.13	7.48	0.05	8.57
	8 July	6: Flowering, anthesis	1016	5.35	10.12	1098.14	2.29	0.85	0.52	0.19	9.31					84.90	3.29	0.05	3.14
	24 July	7: Development of fruit	1300	5.26	7.38	1770.14	11.59	0.86	0.99	0.21	3.36					82.48	6.16	0.05	9.62
	6 August	7: Development of fruit	1534	5.29	16.17	2565.51	16.79	0.87	0.53	0.25	4.87					78.23	7.35	0.04	18.80
	19 August	8: Ripening	1757	4.77	8.05	2910.24	9.11	0.84	0.41	0.22	4.62					66.93	8.37	0.02	51.00
	2 September	8: Ripening	1984	2.14	6.48	2477.75	11.03	0.73	1.29	0.10	12.75					16.15	8.42	−0.21	−4.53
2016	14 June	1: Leaf development	561	0.31	49.49	32.26	25.04	0.37	7.54	−0.14	−14.65	0.29	11.93	−0.20	−11.90	23.00	4.16	−0.15	−3.94
	27 June	3: Stem elongation	739	1.14	41.16	130.56	66.27	0.56	12.39	−0.04	−148.53	0.57	10.75	−0.05	−46.35	32.34	2.30	−0.10	−4.12
	13 July	6: Flowering, anthesis	1011	3.88	10.44	362.86	76.27	0.79	4.07	0.12	25.60	0.82	6.48	0.25	33.29	58.26	16.61	0.03	33.49
	21 July	6: Flowering, anthesis	1143	4.03	13.61	650.56	23.61	0.83	1.64	0.17	8.54	0.87	3.36	0.37	16.33	64.15	12.14	0.04	12.13
	3 August	7: Development of fruit	1368	3.85	13.57	1102.88	2.59	0.84	0.89	0.19	6.39	0.97	0.31	0.64	2.35	66.63	3.59	0.03	15.82
	19 August	7: Development of fruit	1628	4.35	13.22	1533.02	27.67	0.86	1.16	0.25	2.41	0.94	0.23	0.45	3.23	75.88	6.00	0.05	11.96
	26 August	8: Ripening	1745	3.65	6.36	2074.10	25.91	0.81	0.40	0.16	6.39	0.88	1.18	0.35	9.87	72.72	6.45	0.05	17.76
	7 September	8: Ripening	1952	3.64	7.07	2535.06	20.61	0.73	0.07	0.05	7.47	0.80	2.85	0.18	17.95	74.96	4.53	0.05	10.79
	23 September	8: Ripening	2160	1.48	58.53	2924.51	11.04	0.57	0.64	−0.07	−3.51	0.59	5.40	−0.09	−35.13	34.81	19.75	−0.09	−41.99

LAI: leaf area index (m² m⁻² soil); DAGB: dry aboveground biomass (g m⁻²); NDVI_Landsat: normalized difference vegetation index obtained from Landsat; VARI_Landsat: visible atmospherically resistant index from Landsat; NDVI_Sentinel: normalized difference vegetation index obtained from Sentinel; VARI_Sentinel: visible atmospherically resistant index from Sentinel; GCC: green canopy cover (%); VARI_UAV: visible atmospherically resistant index from unmanned aerial vehicle (UAV); $\overline{x}$: mean value; CV: coefficient of variation (%). ^a Bundesanstalt Bundessortenamt and Chemical Industry scale (BBCH) scale [27].

Table 4

Trained models using the 2011, 2012, and 2015 field and remote sensing data.

Variable	Predictors	Equation	R2adj	RMSE
LAI	NDVILandsat, GDD	LAI = −2.4829 + 9.1596NDVILandsat + 2.1015*10^-5GDD	0.723	1.01
	VARILandsat, GDD	LAI = 3.0711 + 11.68VARILandsat−0.00012514GDD	0.613	1.19
	GCCUAV, GDD	LAI = −0.75607 + 5.4966GCCUAV + 0.00095473GDD	0.931	0.468
	VARIUAV, GDD	LAI = 2.8943 + 14.353VARIUAV + 0.0010841GDD	0.853	0.682
DAGB	NDVILandsat, GDD	DAGB = −174.24–1317.3NDVILansat + 2.3997GDD	0.770	589
	VARILandsat, GDD	DAGB = −1001.9–2106.3VARILansat + 2.4677GDD	0.779	578
	GCCUAV, GDD	DAGB = −941.17 + 426.42GCCUAV + 2.0389GDD	0.746	615
	VARIUAV, GDD	DAGB = −685.97 + 785.29792.66VARIUAV + 2.0598GDD	0.741	621

LAI: leaf area index (m² m⁻² soil); NDVI_Landsat: normalized difference vegetation index obtained from Landsat; GDD: growing degree days; VARI_Landsat: visible atmospherically resistant index from Landsat; GCC: green canopy cover (%); VARI_UAV: visible atmospherically resistant index from unmanned aerial vehicle (UAV); DAGB: dry aboveground biomass (g m⁻²); R²_adj: adjusted coefficient of determination; RMSE: root mean squared error (m² leaf m⁻² soil for LAI and g m⁻² for DAGB).

Table 5

Trained and calibrated models using the 2016 field and remote sensing data from Sentinel.

Variable	Predictor	Equation	R2adj	RMSE
LAI	NDVISentinel, GDD	LAI = −1.9427 + 6.5039NDVISentinel	0.871	0.54
	VARISentinel, GDD	LAI = 1.8985 + 4.906VARISentinel	0.754	0.747
DAGB	NDVISentinel, GDD	DAGB = −1111.4 + 1.7402GDD	0.944	244
	VARISentinel, GDD	DAGB = −1111.4 + 1.7402GDD	0.944	244

LAI: leaf area index (m² m⁻² soil); NDVI_Sentinel: normalized difference vegetation index obtained from Sentinel; GDD: growing degree days; VARI_Sentinel: visible atmospherically resistant index from Sentinel; DAGB: dry aboveground biomass (g m⁻²); R²_adj: adjusted coefficient of determination; RMSE: root mean squared error (m² leaf m⁻² soil for LAI and g m⁻² for DAGB).

Table 6

Trained multilinear regression model (MRM) model using the 2016 field and remote sensing data from Sentinel.

Variable	Predictor	Equation	R2adj	RMSE
DAGB	NDVISentinel, GDD	DAGB = −1092.6–29.079NDVISentinel + 1.7415GDD	0.942	248
	VARISentinel, GDD	DAGB = −1111.4 + 15.86VARISentinel + 1.7407GDD	0.942	248

NDVI_Sentinel: normalized difference vegetation index obtained from Sentinel; GDD: growing degree days; VARI_Sentinel: visible atmospherically resistant index from Sentinel; DAGB: dry aboveground biomass (g m⁻²); R²_adj: adjusted coefficient of determination; RMSE: root mean squared error (m² leaf m⁻² soil for LAI and g m⁻² for DAGB).