Analysis of Vegetation Indices to Determine Nitrogen Application and Yield Prediction in Maize (Zea mays L.) from a Standard UAV Service

Abstract

The growing use of commercial unmanned aerial vehicles (UAV) and the need to adjust N fertilization rates in maize (Zea mays L.) currently constitute a key research issue. In this study, different multispectral vegetation indices (green-band and red-band based indices), SPAD and crop height (derived from a multispectral compact camera mounted on a UAV) were analysed to predict grain yield and determine whether an additional sidedress application of N fertilizer was required just before flowering. Seven different inorganic N rates (0, 100, 150, 200, 250, 300, 400 kg·N·ha⁻¹), two different pig slurry manure rates (Ps) (150 or 250 kg·N·ha⁻¹) and four different inorganic-organic N combinations (N100Ps150, N100Ps250, N200Ps150, N200Ps250) were applied to maize experimental plots. The spectral index that best explained final grain yield for the N treatments was the Wide Dynamic Range Vegetation Index (WDRVI). It identified a key threshold above/below 250–300 kg·N·ha⁻¹. WDRVI, NDVI and crop height showed no significant response to extra N application at the economic optimum rate of fertilization (239.8 kg·N·ha⁻¹), for which a grain yield of 16.12 Mg·ha⁻¹ was obtained. This demonstrates their potential as yield predictors at V12 stage. Finally, a ranking of different vegetation indices and crop height is proposed to overcome the uncertainty associated with basing decisions on a single index.

1. Introduction

Nitrogen (N) fertilization of maize (Zea mays L.) is an important research topic. Nitrogen, together with genetic improvement, is one of the most important factors affecting production and can account for up to 30% of the total cost of producing maize [1]. N fertilization is universally accepted as a key input for increasing maize grain yields and optimizing economic returns [2]. In many irrigated Mediterranean areas, maize is one of the most important field crops. In the Ebro valley (NE Spain), grain yields commonly range from 12 to 15 Mg·ha⁻¹, with a total plant nitrogen uptake of 250–300 kg·ha⁻¹ [3]. Farmers usually decide N application rates on the expected crop N uptake which is, in turn, based on yield goals. However, they do not often consider the possible effect of having high N levels in the soil prior to planting, conditions which are common in many crop-growing areas [4]. In irrigated maize, most of the N fertilizer tends to be applied at planting or during the earliest stages of crop growth, as this simplifies crop management. This may constitute a problem since N uptake does not all occur at the same time. Under favourable soil moisture conditions, approximately one-third of the total N uptake occurs after pollination [5]. Consequently, consideration should be given to N applications via sidedress, applied at or near flowering (VT stage or tasseling, [6]), with appropriate doses of N being applied at planting and/or at lay-by (e.g., V7–8, 7–8 leaves with visible leaf collars) in order to maintain yield potential [7].

Over-fertilization can occur as a consequence of N fertilization in early stages of crop growth. This does not increase grain yield but, instead, wastes fertilizer, increases costs, and can cause nitrate pollution [8]. This problem has led farmers, scientist and politicians to explore how to improve N efficiency, reduce N inputs and prevent water and soil pollution associated with maize production [9,10,11].

In recent years, the use of chlorophyll meters (CM), which measure leaf chlorophyll content to estimate N nutrition status, has increased among researchers and farmers [10,11]. Hawkings et al. [12] found that the adjusted R² of the relationship between CM readings and the nitrogen rate difference (ND) for the economic optimum nitrogen rate (EONR) was 0.76 for a maize-maize rotation. However, despite the good correlations, these methods fail to capture the spatial variability that is often present within plots. Data acquisition also tends to be time-consuming and can be hindered by a series of practical limitations [8,9]. One alternative to ground-based measurements involves image-based (satellite or airborne) remote-sensing. This technique, either using active or passive sensors, has been recognized as a potential tool for both spatial and temporal improvement of N management in field crops [13,14], and can also be used to detect N deficiencies in maize [15]. Cilia et al. [16], using multispectral airborne images, created a variable rate N fertilization map based on the difference between actual and optimal crop N content. This map of maize N content also related well to the real maize N content obtained using traditional destructive measurement techniques (R² = 0.70). In this respect, multispectral satellite images, such as Aster and QuickBird, and manned-flight airborne images have been used to assess irrigated maize N status at V12 (12 leaves with visible leaf collars) and later stages of crop growth, as well as to determine field variability for in-season N management in order to complement ground-based measurements [9,10,17]. In fact, most previous studies have shown that both real and false colour images acquired between growth stages V7 and VT could be used to predict N deficiencies and N requirements in maize [18]. Combining CM readings and aerial and/or satellite remote sensing images would therefore seem to offer a practical solution to in-season site-specific N applications in large fields [19].

Image acquisition at stage V12 is generally preferred to other alternatives because the observed maize N uptake at this stage tends to be about 40% of the total [6], and crop response to N fertilizer is high if N deficiencies are detected. In addition, when taking aerial images at stage V12, strong background reflectance from the soil is minimized. This is one of the most challenging obstacles to detecting maize N deficiencies in the early stages of crop growth [18]. At the same time, tassel colour interference, which can occur if the images are taken at later stages, is also avoided. Readings for earlier vegetation stages can usually be discarded since only weak correlations are found (R² ≤ 0.29) for the prediction of optimum N rates [20].

Despite the apparent advantages offered by remote sensing from satellites and manned aircraft, the cost of obtaining high-resolution multispectral images for relatively small areas is considered an important drawback [10,21]. At present, this can be overcome by using unmanned aerial vehicles (UAV), mounted with multispectral cameras. Image acquisition with UAV can be deployed quickly and repeatedly, meaning lower costs, greater flexibility in terms of flying heights and mission timing, and higher spatial resolutions [22]. In recent years, UAV-based research has been carried out to monitor vegetation for agricultural purposes [16]. Most of these applications have been possible due to the miniaturization of multispectral and thermal cameras. However, radiometric and geometric calibrations are required to provide images that are similar to those available from traditional satellite-mounted sensors [23].

Several studies on the use of UAV in the assessment of N status in maize have been published. These have mainly focused on standardizing NDVI (Normalized Difference Vegetation Index) values for their use as part of an N sufficiency index (the NDVI reading divided by the NDVI value of a corresponding well-fertilized N field) [24], and on comparing ground-sensor measurements with hyperspectral images [11]. In [11], indices based on UAV hyperspectral imagery were used to calculate greenness, chlorophyll and photochemical indices. These indices were found to be as reliable as ground-level measurements for assessing crop nitrogen (N) status. This finding was also in line with the work of Isla et al. [17], who investigated the N nutritional status of maize using multispectral data acquired from an aircraft. Scharf and Lory [14] also showed that maize colour measured using aerial photography could be used to predict N sidedress requirements. Correlations between colour and the EONR ranged from 0.60 to 0.79 after the removal of soil pixels.

Other multispectral indices have been reported as being particularly useful for assessing maize N status. McMurtrey et al. [25] reported that as N deficiencies increased, leaf reflectance increased in the green band (0.55 µm), decreased in the NIR (Near-Infrared) (0.70 µm) and remained almost unvaried in the red band (0.67 µm). Along the same lines, Bausch and Duke [8] proposed an N reflectance index for maize based on the Green Ratio Vegetation Index (GRVI), which correlated highly with the N sufficiency index (average SPAD reading for a given treatment divided by the average SPAD reading for a well-fertilized N field; with SPAD being a device for indirect chlorophyll measurement) after stage V11 (11 leaves with visible leaf collars). More recently, other authors have reported the greater sensitivity of the GRVI than red-based indices for assessing N deficiencies in maize [26]. Other green-based vegetation indices have also been qualified as particularly useful for assessing maize N status at V12 or later growth stages [9]. For example, Schlemmer et al. [27] and Li et al. [28] showed the significance of the green and red-edge bands for estimating chlorophyll and N content in maize. The GNDVI (Green Normalized Difference Vegetation Index) has also been considered useful for assessing leaf chlorophyll variability when the leaf area index is moderately high [29]. In other cases, green-based indices have shown high correlations with maize grain yield, explaining 86% of the observed variance [17].

In view of the growing application of UAV services and of the importance of N fertilization in maize production, the purpose of the present study was to analyse the potential utility of different multispectral vegetation indices in conjunction with crop height to support decisions regarding the need to apply N fertilizer just before flowering (V12). The study involved the use of a standard UAV service for acquiring multispectral images with a broad-band compact camera, from which crop height data were automatically derived. It explores the possibilities of using this type of service in commercial farms on a day-to-day basis.

2. Materials and Methods

2.1. Study Area

The experimental field was located at the IRTA Research Station in Gimenells (Lleida, Northeast Spain, 41°65′N, 0°39′E) and had an area of 110 × 130 m² (Figure 1). Soils were well-drained, had no salinity problems, and were characterized by the presence of a petrocalcic horizon at a depth of 80–100 cm, being classified as Petrocalcic Calcixerept [30]. The area has a semi-arid climate with a mean temperature of 19.1 °C and low precipitation during the maize growing season (192 mm) [31]. Irrigation is therefore required to achieve high grain yields. The field was irrigated using a sprinkler irrigation system, which provided approximately 750 mm of water (with no appreciable nitrate content) over the maize growing season. Conventional tillage was applied, which included disc ploughing and cultivation to a depth of 25–30 cm. The study field was divided into 45 experimental plots, each of 10 × 15 m². These plots were sown with maize of the variety PR33Y72 (FAO cycle 600) on 10 April 2014.

A pre-emergence herbicide (S-Metolachlor 40% and Terbutilazine 18.75%) was applied at 3 L·ha⁻¹ to control weeds. At post-emergence, 1 L·ha⁻¹ of Dimethylamine salt of dicamba 48.2% (3.6-dichloro-o-anisic acid) and 0.75 L·ha⁻¹ Nicosulfuron 6% were applied to control Abutilon theophrasti M. and Sorghum halepense L., respectively.

The maize was harvested on 29 October 2014. Grain yield was determined by harvesting two complete central rows (1.5 × 10 m²) using a small-plot combine. A grain sample of 250 g was taken from each plot to determine moisture content (GAC II, Dickey-John, Auburn, IL, USA) and adjust grain yield to 14% moisture.

It should be noted that the plots had been continually used throughout the previous 12-year period as an experimental field for studies related to maize N fertilization. It should also be noted, in this respect, that the same treatments were applied to each individual plot as in the previous years. As a result, the effect of the different N treatments would have been expected to be more evident.

2.2. Fertilizer Treatments

Seven different inorganic N rates were applied: 0, 100, 150, 200, 250, 300 and 400 kg·N·ha⁻¹. Two different pig slurry manure rates (Ps): 150 and 250 kg of organic N·ha⁻¹ were used in the study (Figure 1) and four different inorganic-organic combinations were also applied: N100Ps150, N100Ps250, N200Ps150, and N200Ps250.

Each N treatment included four replications (except for the treatments with Ps250, with three replications) under a split plot design. Pig slurry (Ps) applications were done before sowing and the slurry was then ploughed into the soil between 3–5 h after application to reduce ammonia (NH₃) volatilization losses. Inorganic N fertilization (33.5% N, as ammonium nitrate) was applied by hand, in two equal parts: 50% at the first sidedress (stage V3–V4: 3–4 leaves with visible leaf collars) and 50% at the second sidedress (stage V6–V7: 6–7 leaves with visible leaf collars) [31]. In the case of the N400 treatment, a third sidedress was applied at stage V12 (when the UAV images were acquired). In this case, N application was distributed differently (37.5% at stage V3–V4, 37.5% at V6–V7 and 25% at V12) with a view to reducing the risk of pollution by nitrate leaching, which is associated with high rates of N application. Additionally, phosphorus and potassium fertilizations were applied before planting, at rates of 150 kg·P₂O₅·ha⁻¹ and 250 kg·K₂O·ha⁻¹, to ensure no deficit of either of these elements. This type of fertilization had been performed each 2 years based on previous soil analysis.

2.3. Remote and Proximal Sensing Data Acquisition and Analysis

An aerial survey was carried out with the Atmos-6 UAV (CATUAV, Moià, Catalonia, Spain) (Figure 2). This drone has a wingspan of 1.80 m, a length of 1.29 m and a payload of 500 g. The survey was conducted on 30 June 2014, at 10 h (GMT), with maize at V12 stage. The flight height was 180 m above ground, with a speed of 38 km h⁻¹. The time of flight was very short since the area to capture was only about 0.8 ha, and irradiance conditions did not vary during image acquisition. The images were acquired using a VEGCAM-Pro camera, with a 14 Mp Foveon X3 image sensor. This camera has a total weight of 307 g and a size of 110 × 85 × 78 mm. It works in three wide spectral bands: green (525–575 nm), red (615–685 nm) and near infrared (755–805 nm), with a radiometric resolution of 8 bits/pixel (with a pixel value range of 0–255). Other characteristics of the camera include: objective 16.6 mm/F4, sensor size 20.7 × 13.8 mm², effective pixels 2650 × 1768 × 3 layers, instantaneous field of view (IFOV) 0.0265° × 0.0265° deg.

The photos were acquired with a horizontal overlap of at least 60% to allow the use of stereoscopy to compute the elevation in each pixel. The images were then rectified and mosaicked with the aid of Pix4D software (Pix4D SA, Lausanne, Switzerland) to produce images at a spatial resolution of 0.15 m. The geometry of the camera was accurately calibrated using a calibration panel and the RapidCal software (PIEneering, Helsinki, Finland), yielding standard errors for the principal point of ±0.0122 mm and ±0.0098 mm in X and Y respectively. The photos were georeferenced on the basis of ground control points selected on a 1:2500 scale orthophoto, produced by the Cartographic and Geologic Institute of Catalonia. Due to the small size of the experimental field (0.8 ha), only three sequential photos were needed to produce the mosaic, and artifacts were not produced. Reflectance values were computed for each band by dividing the pixel values by those of a calibrated diffuse Spectralon reflectance target (Labsphere, North Sutton, NH, USA).

Three vegetation indices were computed based on the reflectance values: NDVI (Normalized Difference Vegetation Index, Equation (1), [32]), GRVI (Green Ratio Vegetation Index, Equation (2), [20]), and WDRVI (Wide Dynamic Range Vegetation Index, Equation (3), [33]).

$$NDVI=\frac{NIR-Red}{NIR+Red}$$

(1)

$$GRVI=\frac{NIR}{Green}$$

(2)

$$WDRVI=\frac{\left(\alpha ·NIR-Red\right)}{\left(\alpha ·NIR+Red\right)}$$

(3)

where NIR is the reflectance of the near infrared light, Red is the reflectance of the red light, Green is the reflectance of the green light and α a weighting coefficient that can vary from 0.1 to 0.2. The WDRVI was created to increase correlations with the vegetation fraction for crops such as wheat, soybean and maize, thereby enabling a more robust characterization of the physiological and phenological characteristics of the crop [33]. In the present study, we used α = 0.1 because of its better fit to N dose and maize yield.

In addition, the crop height in each pixel was calculated based on an intensive photogrammetric analysis of stereopairs with Pix4D software. This made possible the creation of a digital surface model (DSM) of the experimental plots. A digital ground model (DGM) based on these data was also built by selecting bare soil pixels around the experimental plots and capturing their height from the DSM. By connecting the different bare soil pixels, it was then possible to create a triangulated irregular network representing ground elevation. Finally, maize plant height was calculated by subtracting DGM from DSM.

In-field crop height measurements were taken at the VT stage (tasseling) (10 to 15 days after the UAV flight). These measurements were then used to calculate their correlation with those obtained via the photogrammetric process. For this, the heights of 5 plants were measured in each experimental plot using a tape measure. These height measurements were repeated until tassel insertion, in order to avoid the possibility of differences in height caused by the length of the tassel. According to Duncan et al. [34], tassel size varies with both plant population and variety. Linear regression analysis was performed between the image and in-field height, yielding an R² coefficient of 0.82 (RMSE = 0.15 m and p-value < 0.001). This indicates a good correlation between the two types of measurement and that the height values obtained from the digital height model could be used as a good estimator of crop height at the moment of image acquisition.

Non-destructive chlorophyll readings were also taken from plant leaves in the experimental plots at VT stage, with a view to comparing the measurements with the spectral indices and correlate them with yield. These measurements were taken using a small, lightweight, portable, hand-held meter (SPAD-502 indirect chlorophyll meter; Minolta Corp, Ramsey, NJ, USA). SPAD values calculate relative chlorophyll content based on the amount of light transmitted by the leaves at two different wavelengths: red (650 nm) and near infrared (940 nm). In agriculture, the SPAD meter is often used to improve N management and increase yields by predicting N status and determining fertilization requirements [35]. In this case, 5 plants were sampled from each plot and three readings were taken from each selected ear leaf: from the base, middle and top of each plant. The 15 measurements taken from each experimental plot were then averaged to obtain the plot SPAD value.

2.4. Statistical Analysis

To analyse the spectral indices and crop height according to the different N fertilizer treatments, fifty points were randomly sampled within each individual plot, excluding a 1 m buffer from the borders. In addition, the points that lay on bare soil were moved to a nearby plant to avoid measurements in areas without crop. A multiple comparison analysis using Tukey’s Honest Significant Difference (HSD) test at a significance level of 0.05 was then carried out. For each N treatment, this analysis compared the mean values of the spectral indices (included SPAD) and crop height, distinguishing different significant groups. JMP Pro 12 (SAS Institute, Cary, NC, USA) statistical package was used for the statistical analysis.

Linear-plateau models were fitted for spectral indices, crop height and yield with the N fertilization treatments. In this way, the economic optimum rates of fertilization, as well as the saturation point for N fertilization of spectral indices and crop height were identified by locating the intersection of the two lines [36].

Linear and quadratic regression analyses were also carried out between the mean values for the spectral indices and crop height from the different treatments and yield. These were used to determine the best variable to estimate yield at the V12 stage.

3. Results and Discussion

3.1. Multiple-Rank Analysis of Spectral Indices, Crop Height and Yield

Table 1. Mean values of the vegetation indices and crop height (n = 50 for each treatment replication), SPAD (n = 15 for each treatment replication) and yield (n = 4) for the different nitrogen (N) treatments applied on the experimental plots. Tukey’s HSD test: different letters indicate homogeneous groups with respect to the mean differences at a p-value of <0.05. The total amount of N (kg·ha⁻¹) applied up to vegetation stage V12 is given in brackets.

N Treatment	NDVI	WDRVI	GRVI	SPAD	Crop Height (m)	Yield (Mg·ha−1)
N0	0.451 d	−0.574 f	1.253 f	30.06 b	0.81 i	3.16 e
N100 (100)	0.942 c	0.619 e	1.355 e	48.68 a	1.65 h	9.79 cd
N150 (150)	0.963 b	0.747 d	1.366 cde	51.76 a	1.85 fg	12.02 bcd
N0Ps150 (150)	0.926 c	0.584 e	1.362 cde	46.30 a	1.88 fg	8.70 de
N200 (200)	0.974 ab	0.809 bcd	1.367 cde	56.30 a	1.88 fg	15.04 abc
N250 (250)	0.967 ab	0.781 cd	1.363 de	56.71 a	1.86 fg	16.17 ab
N0Ps250 (250)	0.971 ab	0.808 bcd	1.373 abcd	54.26 a	1.99 cde	15.14 abc
N100Ps150 (250)	0.975 ab	0.834 abcd	1.377 abcd	56.30 a	1.99 cd	14.65 abc
N300 (300)	0.971 ab	0.793 cd	1.366 cde	56.23 a	1.82 g	15.40 ab
N400 (300) *	0.979 ab	0.846 abc	1.367 bcde	**	1.93 def	17.64 a
N100Ps250 (350)	0.986 a	0.894 ab	1.385 a	56.48 a	2.11 ab	17.00 ab
N200Ps150 (350)	0.987 a	0.900 a	1.383 ab	56.89 a	2.03 bc	17.56 a
N200Ps250 (450)	0.987 a	0.899 ab	1.382 abc	57.36 a	2.16 a	16.98 ab

* A total of 300 kg·N·ha⁻¹ was applied up to stage V12 (image acquisition). ** Not measured.

shows the mean values of the vegetation indices (n = 50 for each treatment replication), crop height (n = 50), SPAD (n = 15 for each treatment replication) and yield (n = 4) for the different N treatments applied on the experimental plots. The table also presents the mean values of the variables analysed. The treatments are sorted according to the total amount of N applied until stage V12. Treatments not connected by the same letter are significantly different at a p-value of <0.05. Figure 3 provides a visual comparison of the vegetation indices and crop heights in the experimental fields.

For all the considered variables, the N fertilizer treatments were grouped in three to nine homogeneous groups. Of these, N0 (the control treatment) presented the lowest values for all of the variables and also for grain yield (3.16 Mg·ha⁻¹). The N0 treatment was clearly identified as being totally different from the other N treatments, with no nitrogen applications in these plots over the 12 continuous years of the experiment. The treatments involving applications between 0 and 150 kg·N·ha⁻¹ were separated into several different homogeneous groups, which in most cases exhibited clear N deficiencies in their spectral indices. The N100 and N0Ps150 treatments were associated with the lowest grain yields and the lowest index values after the N0 treatment. The N150 treatment resulted in higher grain yield than the N0Ps150 treatment due to the greater efficiency of applying N at sidedress rather than at planting [37,38].

The treatments which combined inorganic and organic N, with a total amount of applied N of between 350 and 450 Kg·N·ha⁻¹ (N100Ps250, N200Ps150 and N200Ps250), were at the top of the ranking. This seems to indicate that up to V12 stage, the plots fertilized with a combination of organic and inorganic N performed better in terms of plant vigour than those only fertilized with inorganic N. The reason for this may be a progressive mineralization of organic N, which would have facilitated N availability during all development stages [39,40]. The experimental plots in which these treatments were applied also produced the highest grain yields, although not exclusively, as the N200, N250, N0Ps250, N100Ps150, N300 and N400 were also classified in homogeneous group “a”, the one with the highest yields.

The results obtained for the N400 treatment require particular attention. Based on the total amount of N applied, this treatment might have been expected to be in the upper part of the homogeneous group classification (group “a” or “ab”) for the different variables considered. However, it was not, even though the N400 treatment did produce the highest grain yield (17.64 Mg·ha⁻¹). The explanation for this lack of correspondence is related to the total amount of N applied on the experimental plots before image acquisition. Sidedress applications were completed at stage V6-V7 for all treatments except the N400, with the final sidedress for the N400 plots being administered at stage V12. As result, only 75% of the N fertilizer (300 kg·N·ha⁻¹) was applied in those plots before the image acquisition date.

According to the vegetation indices and crop height (Table 1), a total application of 250–300 kg·N·ha⁻¹ corresponds to the threshold value. Above this quantity, the mean values of the analysed variables were significantly higher than for the other N treatments (except for the SPAD, which showed a homogeneous response to the different treatments). This, together with the final performance of the N400 treatment, seems to indicate the possibility at V12 stage of determining whether a supplementary N sidedress application is required to achieve a higher yield.

Among the variables analysed (Table 1), the GRVI performed better than the NDVI at identifying the treatments equal to or below 300 kg·N·ha⁻¹. The latter showed an increase in values that was in line with the total amount of applied N. However, this response was either more homogeneous or greater than for applications of 200 kg·N·ha⁻¹,yet without showing a clear discrimination between the treatments occupying the highest positions in the ranking. Nevertheless, GRVI did not clearly differentiate between the N150 and N300 treatments, and this created a certain degree of ambiguity. In addition, GRVI values for the N250 and N300 treatments did not directly correlate with the total amount of applied N. These results were in line with those published in other research work. Isla et al. [17] reported a better performance of green-based vegetation indices in determining maize vigour due to problems of saturation associated with NDVI for some types of vegetation during their later stages of growth. For example, indices such as the GNDVI have been considered more useful for assessing leaf chlorophyll variability when the leaf area index is moderately high [29]. Xiang and Tian [22] also reported that the GNDVI and GRVI offered the best ways to identify three different N treatments over the whole maize growing season, with the greatest differences in index values being observed during the V6–V8 stages.

The WDRVI was the spectral index that best distinguished between N applications above 300 kg·N·ha⁻¹, but also displayed a certain degree of ambiguity when it came to distinguishing between treatments within the range between N200 and N300. Only one exception to this general rule was observed; this occurred with the N100Ps150 treatment. This was included in the group of treatments classified as “a”, although its results also overlapped with those of other groups. According to this finding, mean WDRVI values for maize fields of less than 0.89 at V12 would be associated with an improvement in grain yield following the application of a third N sidedress around the VT stage.

In the case of crop height measured from the UAV photogrammetric survey, it was observed that the combined inorganic and organic N treatments at doses of 350 and 450 kg·N·ha⁻¹ were clearly associated with the largest crop heights. Nine significant different responses were identified among the 13 N treatments. The treatments with total N rates of between 150 and 300 kg·N·ha⁻¹ were not classified in ascending order. Nevertheless, a clear differentiation was achieved between treatments that included organic N and those that only contained inorganic N. This was probably due to the effect of higher N availability in the organic plots during the earlier vegetative stages, prior to the sidedress application, which resulted in faster maize plant growth.

With respect to grain yield, the fertilizer treatments N100Ps250, N200Ps150 and N400 showed the greatest response to N (with grain yields equal to or above 17 Mg·ha⁻¹), and were classified in groups “ab” (N100Ps250) and “a” (N200Ps150 and N400). The N treatments with total N content of between 200 and 300 kg·ha⁻¹ and the treatment with the highest N application rate (N200Ps250) were classified as either “ab” or “abc” and produced grain yields of between 15 and 17 Mg·ha⁻¹.

To clarify the type of relationships between the N rates and the analysed variables, linear-plateau models were fitted (Figure 4). In this way, it was possible to identify the economic optimum rates of fertilization, as well as the saturation point for N fertilization of vegetation indices and crop height at V12, SPAD at VT and grain yield.

The predicted economic optimum rate of fertilization was determined at 239.8 kg·N·ha⁻¹ (Figure 4F). Except for the GRVI, the other vegetation indices and crop height showed no response to higher N rates above the economic optimum rate. Particularly useful to predict yield response in relation to N rates at V12 stage were the NDVI, WDRVI and crop height, with no response to N rates higher than 247.5, 243.1 and 243.0 kg·N·ha⁻¹, respectively. For SPAD, saturation was determined at a lower N rate (203.8 kg·N·ha⁻¹). This could imply an underestimation of the N requirements of maize, which could be translated into a reduction in final grain yield. The opposite was observed for the GRVI, where N requirements may well be overestimated, with an increase in GRVI values at higher N rates (up to 362.2 kg·N·ha⁻¹). This was probably due to an increase in greenness as observed through the index, which was not translated into final grain yield (Figure 4B).

As for grain yield (Figure 4F), there would be no interest in increasing N application above 239.8 kg·N·ha⁻¹, partly because of potential environmental problems but also because of the subsequent reduction in the economic margin. However, higher N rates did present higher yields, suggesting the potential to increase grain yield, since yield values for N treatments N100Ps250, N200Ps150 N200Ps250 were above the linear-plateau which was reached at 16.12 Mg·ha⁻¹. This suggests the potential possibility of increasing yield by increasing N fertilization, something which would be of particular interest in maize-growing areas like the Ebro valley which have the appropriate environmental conditions and irrigation facilities to attain such high yields.

Table 2. Nitrogen fertilization treatments ordered by the mean values for each vegetation index and crop height at stage V12. The values are ordered from 1 to 13 (1 = lowest vegetation index and 13 = highest vegetation index). The treatments are ordered according to the “Rank” column, which orders the sum of the ranks. The total amount of N (kg·ha⁻¹) applied up to vegetation stage V12 is given in brackets.

N Treatment	NDVI	GRVI	WDRVI	Crop Height	Sum	Rank
N0 (0)	1	1	1	1	4	1
N100 (100)	3	2	3	2	10	2
N0Ps150 (150)	2	3	2	6	13	3
N150 (150)	4	5	4	4	17	4
N250 (250)	5	4	5	5	19	5
N300 (300)	7	6	6	3	22	6
N200 (200)	8	8	8	7	31	7
N0Ps250 (250)	6	9	7	9	31	8
N400 (300)	10	7	10	8	35	9
N100Ps150 (250)	9	10	9	10	38	10
N100Ps250 (350)	11	13	11	12	47	11
N200Ps150 (350)	12	12	13	11	48	12
N200Ps250 (450)	13	11	12	13	49	13

was created to offer a summary of the results presented above, and shows the treatments ordered according to the mean values of the indices and crop height. These values are ranked (from lowest to highest mean value) for the different N treatments. For each variable, the lowest rank is equal to 1 and the highest rank is equal to 13. The “Sum” column in Table 2 indicates the sum of the ranks in each row (or for each N treatment). In this case, the sum of ranks represents an order of magnitude of each N treatment with respect to its spectral response (vegetation index) and crop height. It is useful to compare the N treatments not only according to the values assigned to each vegetation index or crop height, but considering them all together. The “Rank” column shows the final ranking according to the Sum from the lowest to the highest value.

The ranks shown in Table 2 indicate that the treatments with the highest responses in the spectral indices and in terms of crop height were those that combined inorganic and organic N, with total application rates of 350 or 450 kg·N·ha⁻¹ (N200Ps150, N100Ps250 and N200Ps250). These were followed, in descending order, by the treatments containing either only inorganic N or only organic N; the exceptions were treatments N0Ps250 and N100Ps150, which respectively occupied the 8th and 10th positions in the ranking. The control treatment (N0) occupied the lowest position.

Two treatments providing N from pig slurry manure (N0Ps250 and N100Ps150) were ranked in high positions, but did not finally achieve the corresponding yield rank (7th and 5th respectively). This could have been due to differences in the timings of the N applications. The organic applications were applied before the maize was sown and the inorganic applications at sidedress. In addition, there were differences in the amount of N applied with respect to that applied in other treatments at the same time. These pig slurry applications would have conferred optimum conditions for crop development until V12, the stage at which the images were acquired. However, the reduced amount of total N applied (250 kg·N·ha⁻¹) and/or the lack of a third sidedress had a determining influence on the final yields associated with these treatments.

3.2. Relationship between the SPECTRAL Indices and Crop Height with Yield

Table 3. Linear and quadratic regression adjustments between the spectral indices and crop height with yield, considering all the total amounts of N applied in all the different treatments, except N0 and N400.

Vegetation Index	Linear Regression with Yield			Quadratic Regression with Yield
	R2	RMSE Mg·ha−1	p-Value	R2	RMSE Mg·ha−1	p-Value
SPAD	0.88	0.82	<0.001	0.89	1.10	0.0002
NDVI	0.90	0.98	<0.001	0.92	0.92	<0.001
GRVI	0.64	1.86	0.003	0.71	1.78	0.0075
WDRVI	0.92	0.87	<0.001	0.92	0.93	<0.001
Crop height	0.60	1.97	<0.001	0.64	1.98	0.0174

presents the results of linear and quadratic regression analysis comparing the vegetation indices, SPAD and crop height with grain yield. In general, quadratic regression models did not significantly improve yield prediction with respect to the linear models, at the same time adding complexity. The best correlation was obtained with WDRVI (R² = 0.92 and RMSE = 0.87 Mg·ha⁻¹), followed by NDVI and SPAD (R² = 0.90 and 0.88, respectively). It should be noted that SPAD was measured at VT and did not show significant differences with respect to the different N treatments (except for the N0 treatment, which was statistically different from the other N treatments, Table 1).

The GRVI performed worse than the red-based indices (R² = 0.64 and RMSE = 1.86 Mg·ha⁻¹). These results differed from those obtained in other studies, in which green-based indices at V15 (e.g., GNDVI) were highly correlated with maize yield and explained 86% of the variance [17]. Other authors have also pointed to green-based vegetation indices being particularly useful for assessing N status at stage V12 or at later stages of maize growth [9]. Crop height at V12 was not a good indicator of grain yield (R² = 0.60 and RMSE = 1.97 Mg·ha⁻¹). These results partially agree with those of Yin et al. [41], who presented crop height as a good predictor of yield in irrigated maize, varying between R² = 0.52 and 0.86, depending on the year in a three-year experiment.

The WDRVI was the best spectral index obtained from the VEGCAM-Pro camera mounted on the drone in estimating yield at the V12 stage of maize development. This was in line with the outcomes of multiple rank analysis, which identified WDRVI as the best index for discriminating between total N applications at applications above 250–300 kg·ha⁻¹. The results agree with those from other studies, although those were conducted with data at a very different spatial resolution using the Moderate Resolution Imaging Spectroradiometer (MODIS, 250 m pixel resolution). For example, Sakamoto et al. developed a practical method for near real-time prediction of U.S. maize yield based on the WDRVI taken 7–10 days before the corn silking stage [42,43]. These authors found a strong linear correlation with maize grain yield at both field and regional scales. Wang et al. derived phenology-adjusted spectral indices from MODIS data to be used in developing linear regression models with maize yield data [44]. In this case, the peak correlation between the WDRVI and yield was detected 85 days after green-up date (R² = 0.506). The correlation was generally low for NDVI (R² = 0.385) and no obvious peak correlation existed. In other cases, the WDRVI also performed better than the NDVI with MODIS based data, and has been shown to be useful for assessing early stages of plant stress in maize and soybeans [45]. As far as we are aware, the WDRVI has been mainly applied in maize studies using MODIS data with 250 m pixel resolution. In the present study, it has been shown that good results in estimation of N nutritional state and grain yield can also be obtained at a very high resolution with data from UAV-mounted cameras.

4. Conclusions

This study was motivated by the growing demand for commercial UAV services in agriculture and the need to adjust N fertilization rates applied to maize crops. With this in mind, one green-based (GRVI) and two red-based (NDVI and WDRVI) vegetation indices, as well as SPAD and crop height (derived from a photogrammetric process), were analysed in order to determine crop status at stage V12 (12 leaves with visible leaf collars, just before flowering), associated with different amounts of supplied N.

The results obtained led us to conclude that the spectral index derived from the VEGCAM-Pro camera that explained the greatest variability between treatments was the WDRVI. This index had previously only been applied in maize studies at moderate resolution (250 m per pixel) with MODIS data on phenological characterization and yield prediction. In the present case study, at very high spatial resolutions, WDRVI was the best index for distinguishing between treatments with applications above or below 250–300 kg·N·ha⁻¹ and at grain yield prediction at the V12 stage. NDVI and crop height also showed no significant response to extra N application at the economic optimum rate of fertilization, demonstrating their potential as yield predictors at V12 stage. However, SPAD and GRVI either underestimated or overestimated the optimum N rate.

Although there would theoretically be little interest in increasing N application above 239.8 kg·N·ha⁻¹, the study does show a tendency for increased grain yield with higher N rates. This could be of particular interest in maize-growing areas such as the Ebro valley which have the appropriate environmental conditions and irrigation facilities to attain such high yields.

The ranking of the spectral indices and crop height revealed that the treatments that conferred the greatest responses to N fertilization were those which combined inorganic and organic N with applications of 350 or 450 kg·N·ha⁻¹. The proposed ranking system, which is based on the response of the crop to different multispectral indices and crop height, may help to overcome the uncertainty of decision-making based on a single index, such as the NDVI, which is the approach most frequently used in association with data obtained from multispectral cameras.