**Using UAV Collected RGB and Multispectral Images to Evaluate Winter Wheat Performance across a Site Characterized by Century-Old Biochar Patches in Belgium**

#### **Ramin Heidarian Dehkordi 1,\*, Victor Burgeon <sup>1</sup> , Julien Fouche <sup>2</sup> , Edmundo Placencia Gomez <sup>3</sup> , Jean-Thomas Cornelis <sup>1</sup> , Frederic Nguyen <sup>3</sup> , Antoine Denis <sup>4</sup> and Jeroen Meersmans <sup>1</sup>**


Received: 7 June 2020; Accepted: 3 August 2020; Published: 4 August 2020

**Abstract:** Remote sensing data play a crucial role in monitoring crop dynamics in the context of precision agriculture by characterizing the spatial and temporal variability of crop traits. At present there is special interest in assessing the long-term impacts of biochar in agro-ecosystems. Despite the growing body of literature on monitoring the potential biochar effects on harvested crop yield and aboveground productivity, studies focusing on the detailed crop performance as a consequence of long-term biochar enrichment are still lacking. The primary objective of this research was to evaluate crop performance based on high-resolution unmanned aerial vehicle (UAV) imagery considering both crop growth and health through RGB and multispectral analysis, respectively. More specifically, this approach allowed monitoring of century-old biochar impacts on winter wheat crop performance. Seven Red-Green-Blue (RGB) and six multispectral flights were executed over 11 century-old biochar patches of a cultivated field. UAV-based RGB imagery exhibited a significant positive impact of century-old biochar on the evolution of winter wheat canopy cover (*p*-value = 0.00007). Multispectral optimized soil adjusted vegetation index indicated a better crop development over the century-old biochar plots at the beginning of the season (*p*-values < 0.01), while there was no impact towards the end of the season. Plant height, derived from the RGB imagery, was slightly higher for century-old biochar plots. Crop health maps were computed based on principal component analysis and k-means clustering. To our knowledge, this is the first attempt to quantify century-old biochar effects on crop performance during the entire growing period using remotely sensed data. Ground-based measurements illustrated a significant positive impact of century-old biochar on crop growth stages (*p*-value of 0.01265), whereas the harvested crop yield was not affected. Multispectral simplified canopy chlorophyll content index and normalized difference red edge index were found to be good linear estimators of harvested crop yield (*p*-value(Kendall) of 0.001 and 0.0008, respectively). The present research highlights that other factors (e.g., inherent pedological variations) are of higher importance than the presence of century-old biochar in determining crop health and yield variability.

**Keywords:** UAV remote sensing; crop monitoring; RGB imagery; multispectral imagery; century-old biochar

#### **1. Introduction**

Biochar is a carbon-rich material obtained through the pyrolysis of feedstocks to intentionally amend soil [1,2]. In the context of climate-smart agriculture, special attention has been given to biochar application owing to its supposed ability to improve soil fertility and long-term organic carbon storage [3–6]. To date, short-term biochar effects have been well addressed by studying the soil nutrient availability [7,8], soil water holding capacity [9–11], and plant available water content [12,13]. The aforementioned improvements of soil quality could consequently lead to an impact on the crop productivity [14]. Moreover, biochar was found to increase nutrient uptake and crop yield while being incorporated into soil matrix together with fertilizers [15]. The short-term improvement of crop productivity can be attributed to the addition of nutrients related to the biochar application [7]. However, long-term biochar effects on soil and agro-ecosystems properties are still poorly understood. Several authors have reported a higher nutrient availability in the century-old biochar sites originating from earth kiln sites [16,17]. The higher nutrient availability in long-term biochar-enriched soils can be explained by cation exchange capacity that increases over time [18]. Moreover, century-old biochar significantly affects physico-chemical properties of soils in relic hearths of agricultural fields [16].

Across the globe, many efforts have been devoted to unravel biochar effects on aboveground biomass productivity. In many of these studies, a short-term effect of biochar on increasing crop productivity was reported [7,19,20], however, Güereña et al. [21] and van Zwieten et al. [22] showed a reductive influence. Similarly, Hernandez-Soriano et al. [23] reported a long-term biochar impact on increasing aboveground productivity, whereas Mastrolonardo et al. [24] observed a decreased tree ring width in relic charcoal hearths of Belgium. Oak growth depression was also found in aged charcoal hearth patches in the United States [25].

The studies listed above all explored the biochar effects on crop productivity. However, it is highly likely that biochar could affect crop development only during the green-up phase, while the harvested crop yield could remain unaffected [26]. More specifically, long-term biochar effects on crop performance and its development over the entire growing season have received much less attention than biochar effects on total crop productivity.

In recent years, remote sensing has been promoted as a means in precision agriculture for detailed spatio-temporal monitoring of crop dynamics [27,28]. A particular focus at present is the use of satellite imagery for precision agricultural practices due to the images' fine spatial, spectral, and temporal resolution [28]. Moreover, high-resolution unmanned aerial vehicle (UAV)-based imagery can be used to address the data gap in satellite products suffering from cloud coverage [29]. Recent advances in Red-Green-Blue (RGB) and multispectral UAV remote sensing have improved the ability to accurately monitor crop performance and vegetation characteristics at the canopy level, such as canopy cover, leaf area index, and chlorophyll content, allowing for site-specific decision making in the context of precision agriculture [29–32]. In addition, due to remote UAV flights, the experimental surface is not damaged [33], thus enabling a complete crop monitoring over the entire season. Therefore, the spatio-temporal information provided by the UAV images can reveal variability in crop performance due to the presence of century-old biochar.

Every year, agricultural products suffer significant economic losses caused by multiple pests and diseases, such as yellow rust or aphid, which may lead to a shortage of food production globally [34,35] and particularly in Europe [36]. Cameras onboard UAVs are able to illustrate changes in crop dynamics as a function of a particular disease at early onset [37]. The pioneering work of Franke et al. [38] reported that pathogens reduces plant chlorophyll content in the visible and red-edge spectral bands due to necrotic and chlogotic lesions. Su et al. [37] highlighted that the red band of the spectrum has strong potential in discriminating healthy from diseased winter wheat plants. Su et al. [37] also showed that the multispectral optimized soil adjusted vegetation index (OSAVI) has a strong discriminating capability between healthy and diseased winter wheat plants. In addition, UAV imagery is able to detect changes in crop traits linked to disease incidence [39]. Therefore, high-resolution UAV remote sensing can be successfully applied for spatial explicit assessment of crop health as a consequence of the presence of century-old biochar within agricultural soils.

The primary objective of this paper is to evaluate crop performance by studying crop growth and health through the usage of high-resolution RGB and multispectral information, respectively, derived from two UAV platforms. This approach also allowed us to investigate, as a secondary objective, the impacts of century-old biochar on winter wheat crop performance. The available imagery sources were RGB and multispectral data collected using two UAV platforms. UAV remote sensing data in combination with ground-based measurements of crop traits were used to determine the influence of century-old biochar from kiln sites of an agricultural field on winter wheat crop growth. The present study evaluated, for the first time, the impact of century-old biochar on crop performance at the canopy level using UAV imagery, while paving the way to quantify the long-term effects of biochar on winter wheat crop health.

#### **2. Materials and Methods**

#### *2.1. Site Description and Data Acquisition*

The data were collected in a farm of approximately 13 hectares cultivated with winter wheat located in Isnes, Belgium (Figure 1). Winter wheat (*Triticum aestivum* L.) seeds were sown on 6 December 2018. The region is of temperate climate with dry and hot summers. The soil type is a Luvisol characterized by a silt loam texture [40]. Homogenous agricultural treatments were applied throughout the field. The detailed description is available in Heidarian Dehkordi et al. [26], in which the same experimental configuration was followed by eleven 10 × 10 m century-old biochar plots (red squares in Figure 1) together with the corresponding reference plots (blue squares in Figure 1). Ground control point (GCP) targets, with known centroid coordinates through a real-time kinematic GPS, were deployed throughout the farm (Figure 1) to help with accurate geo-referencing of the data (i.e., 3 cm horizontal accuracy).

Two remote sensing datasets (i.e., RGB and multispectral) were acquired during the 2019 growing season (Table 1) using high-resolution sensors onboard two UAV platforms to monitor the crop status over time (Figure 2). RGB images were collected using a DJI Phantom 4 Pro and the multispectral images were obtained using a MicaSense RedEdge-M and downwelling light sensor, onboard a DJI Matrice 100 platform specifically designed by Heidarian Dehkordi et al. [26], capturing blue (475 nm central wavelength), green (560 nm), red (668 nm), red-edge (717 nm), and near-infrared (840 nm) spectral channels with bandwidths of 20, 20, 10, 10, and 40 nm, respectively. The UAV platforms were programmed to fly at 65 and 60 m above ground altitude (AGL) with an airspeed of 5 and 7 ms−<sup>1</sup> for RGB and multispectral sensors, respectively. The forward and sideway image overlaps were 85% and 75% for RGB, and 75% and 75% for multispectral acquisitions, reaching a ground sampling distance (GSD) of 1.8 and 3.7 cm for RGB and multispectral images, respectively.


**Table 1.** An overview of the data acquisition using the unmanned aerial vehicles within the 2019 growing season. Flight times (start time–end time) are expressed in local time.

**2. Materials and Methods** 

effects of biochar on winter wheat crop health.

*2.1. Site Description and Data Acquisition* 

The primary objective of this paper is to evaluate crop performance by studying crop growth and health through the usage of high-resolution RGB and multispectral information, respectively, derived from two UAV platforms. This approach also allowed us to investigate, as a secondary objective, the impacts of century-old biochar on winter wheat crop performance. The available imagery sources were RGB and multispectral data collected using two UAV platforms. UAV remote sensing data in combination with ground-based measurements of crop traits were used to determine the influence of century-old biochar from kiln sites of an agricultural field on winter wheat crop growth. The present study evaluated, for the first time, the impact of century-old biochar on crop performance at the canopy level using UAV imagery, while paving the way to quantify the long-term

The data were collected in a farm of approximately 13 hectares cultivated with winter wheat located in Isnes, Belgium (Figure 1). Winter wheat (*Triticum aestivum* L.) seeds were sown on 6 December 2018. The region is of temperate climate with dry and hot summers. The soil type is a Luvisol characterized by a silt loam texture [40]. Homogenous agricultural treatments were applied throughout the field. The detailed description is available in Heidarian Dehkordi et al. [26], in which the same experimental configuration was followed by eleven 10 × 10 m century-old biochar plots (red squares in Figure 1) together with the corresponding reference plots (blue squares in Figure 1). Ground control point (GCP) targets, with known centroid coordinates through a real-time kinematic

**Figure 1.** (**a**) A schematic layout of the experimental pairs (reference versus century-old biochar plots) in the winter wheat field and distribution of the ground control points. Background image corresponds to the Red-Green-Blue (RGB) orthomosaic captured by the unmanned aerial vehicle on the 16 April 2019. The right magnified windows present an example of a plot (reference plot 5) for visual RGB (**b**) and multispectral weighted difference vegetation index (**c**) monitored on 16 April 2019. **Figure 1.** (**a**) A schematic layout of the experimental pairs (reference versus century-old biochar plots) in the winter wheat field and distribution of the ground control points. Background image corresponds to the Red-Green-Blue (RGB) orthomosaic captured by the unmanned aerial vehicle on the 16 April 2019. The right magnified windows present an example of a plot (reference plot 5) for visual RGB (**b**) and multispectral weighted difference vegetation index (**c**) monitored on 16 April 2019. (840 nm) spectral channels with bandwidths of 20, 20, 10, 10, and 40 nm, respectively. The UAV platforms were programmed to fly at 65 and 60 m above ground altitude (AGL) with an airspeed of 5 and 7 ms−1 for RGB and multispectral sensors, respectively. The forward and sideway image overlaps were 85% and 75% for RGB, and 75% and 75% for multispectral acquisitions, reaching a ground sampling distance (GSD) of 1.8 and 3.7 cm for RGB and multispectral images, respectively.

**Figure 2.** Methodological flowchart for crop traits retrieval from the unmanned aerial vehicle (UAV) imageries in combination with ground-based data. **Figure 2.** Methodological flowchart for crop traits retrieval from the unmanned aerial vehicle (UAV) imageries in combination with ground-based data.

**Table 1.** An overview of the data acquisition using the unmanned aerial vehicles within the 2019

02/22/2019 10:18–11:41 - 03/20/2019 11:02–12:20 12:46–13:40 03/28/2019 11:33–12:59 13:37–14:32 04/16/2019 12:01–13:31 13:49–14:44 04/29/2019 11:00–12:26 12:44–13:54 05/13/2019 11:03–12:36 13:05–14:18 06/24/2019 11:04–12:26 13:16–14:36 The photogrammetric processing of the acquired UAV images was conducted in Pix4Dmapper photogrammetric software (Pix4D S.A., Lausanne, Switzerland). After the aerial triangulation of the UAV images, the geometric processing was optimized in Pix4Dmapper using the coordinates of the GCPs, and dense point clouds and final orthomosaic images were then created. For the multispectral data, the radiometric calibration was applied using MicaSense calibrated reference panels with

**RGB Multispectral** 

growing season. Flight times (start time–end time) are expressed in local time.

The photogrammetric processing of the acquired UAV images was conducted in Pix4Dmapper photogrammetric software (Pix4D S.A., Lausanne, Switzerland). After the aerial triangulation of the UAV images, the geometric processing was optimized in Pix4Dmapper using the coordinates of the GCPs, and dense point clouds and final orthomosaic images were then created. For the multispectral data, the radiometric calibration was applied using MicaSense calibrated reference panels with known reflectance values (MicaSense, Seattle, WA, USA) taken immediately before and after each flight to convert the raw image digital number (DN) to the calibrated reflectance values. In addition, sensor's sensitivity to light, expressed as ISO, values recorded by the downwelling light sensor were used to aid the radiometric calibration. Both RGB and multispectral datasets were converted to GeoTIFF format with the projected geographic coordinates system WGS 1984 UTM Zone 31N.

In addition to the aforementioned remotely-sensed datasets, ground based measurements of crop traits were performed (Figure 2). Crop growth stages were determined on 28 May 2019 according to the stabiologische bundesanstalt, bundessortenamt and chemical industry (BBCH) scale for winter wheat [41]. For each plot, a regular non-destructive sampling method of 5 plants, i.e., 4 in the corners and 1 in the middle to capture the intra-plot variability, was followed. The BBCH code of each plant was visually determined and the median value of the five plants represented the plot growth stage.

Yield measurements were conducted on 20 July 2019. For each experimental plot, an area of approximately 20 m<sup>2</sup> was harvested to determine the yield.

#### *2.2. UAV-Based Crop Monitoring*

Canopy cover was estimated from the acquired high-resolution RGB images (Figure 2) based on the classification method by Heidarian Dehkordi et al. [26]. This algorithm classifies the pixels in the orthomosaic to vegetation and non-vegetation classes according to the spectral differences between the green band and the blue and red bands as below:

$$\left[ (\text{Green} - \text{Blue}) = \alpha \right] \text{and } \left[ (\text{Green} - \text{Red}) = \beta \right] \tag{1}$$

where Blue, Green, and Red show the digital number (DN) of the corresponding spectral channel of the visible spectrum, α and β are the classification threshold values with α representing the spectral difference between the green and red bands and β revealing the spectral difference between the green and red spectral channels. Similar to the finding of Heidarian Dehkordi et al. [26], a threshold value of greater than 20 for both α and β was also reached in the present study through the visual inspection for the best vegetation classification. Canopy cover estimates were used as a potential discriminative parameter between the reference and century-old biochar plots in each acquisition date [26].

Plant height, as an important phenotyping crop trait, has received no attention in the previous studies exploring the impacts of century-old biochar on crop performance. In this study, plant height was estimated from the RGB images (Figure 3) as below:

$$\text{Height} \left(\text{t}\right) = \text{DSM} \left(\text{t}\right) - \text{DTM} \tag{2}$$

where Height is the computed wheat height, DSM and DTM are the UAV-based digital surface model and digital terrain model, respectively, with index t representing the acquisition date. DSMs were generated from the RGB flights for each acquisition date. The RGB flight on 22 February 2019 (Table 1) was executed when the field was completely bare soil and, hence, the UAV-based DSM, revealing the pure terrain surface topography, was used as the DTM for this study. GCP\_2 (Figure 1) was used as the reference levelling point to calibrate the generated UAV-based DSMs at each acquisition date, bringing the elevation values onto the same scale. The created RGB-based vegetation mask (Equation (1)), with GSD of 1.8 cm, was used to retrieve the DSM from the true canopy pixels for the height computation.

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 6 of 24

**Figure 3.** Methodological workflow to derive plants height from RGB images collected by the unmanned aerial vehicle. DTM is the digital terrain model acquired on 22 February when the field was completely bare soil. DSM represents the digital surface model at each acquisition date. Ground control point 2 was used as a reference levelling point to calibrate the DSM images. **Figure 3.** Methodological workflow to derive plants height from RGB images collected by the unmanned aerial vehicle. DTM is the digital terrain model acquired on 22 February when the field was completely bare soil. DSM represents the digital surface model at each acquisition date. Ground control point 2 was used as a reference levelling point to calibrate the DSM images.

In addition to the RGB derivatives listed above, narrow-band vegetation indices given in Table 2 were computed from the multispectral dataset (Figure 2) as a proxy to crop traits allowing for demonstrating spectral changes related to crop structure [39]. Furthermore, the multispectral information derived from the vegetation indices were analyzed to study the impacts of century-old biochar on winter wheat crop spectral information. The normalized difference vegetation index (NDVI) is a well-known proxy of photosynthetic activity and vegetation greenness. Canopy cover and leaf area index (LAI) were accurately predicted based on the weighted difference vegetation index (WDVI) in several studies [42–44]. The normalized difference red edge index (NDRE) was shown as a good proxy of wheat nitrogen concentration [45]. The optimized soil adjusted vegetation index (OSAVI) was found to be a robust indicator of leaf chlorophyll content [46], and the same conclusion was reached for the chlorophyll vegetation index (CVI). The enhanced vegetation index (EVI) was shown to be sensitive to canopy structural parameters such as LAI while adjusting soil background effects and atmospheric resistance using the blue band of the spectrum [47]. The chlorophyll index red (CI-red) and the simplified canopy chlorophyll content index (s-CCCI) are known alternatives of structural indices which are often used for estimating canopy chlorophyll and nitrogen contents [48,49]. The aforementioned vegetation indices were used for identifying the crop spectral differences between reference and century-old biochar plots throughout the growing season. First, WDVI was In addition to the RGB derivatives listed above, narrow-band vegetation indices given in Table 2 were computed from the multispectral dataset (Figure 2) as a proxy to crop traits allowing for demonstrating spectral changes related to crop structure [39]. Furthermore, the multispectral information derived from the vegetation indices were analyzed to study the impacts of century-old biochar on winter wheat crop spectral information. The normalized difference vegetation index (NDVI) is a well-known proxy of photosynthetic activity and vegetation greenness. Canopy cover and leaf area index (LAI) were accurately predicted based on the weighted difference vegetation index (WDVI) in several studies [42–44]. The normalized difference red edge index (NDRE) was shown as a good proxy of wheat nitrogen concentration [45]. The optimized soil adjusted vegetation index (OSAVI) was found to be a robust indicator of leaf chlorophyll content [46], and the same conclusion was reached for the chlorophyll vegetation index (CVI). The enhanced vegetation index (EVI) was shown to be sensitive to canopy structural parameters such as LAI while adjusting soil background effects and atmospheric resistance using the blue band of the spectrum [47]. The chlorophyll index red (CI-red) and the simplified canopy chlorophyll content index (s-CCCI) are known alternatives of structural indices which are often used for estimating canopy chlorophyll and nitrogen contents [48,49].

used to generate a multispectral-based vegetation mask. A threshold of WDVI higher than 0.30 was used to retrieve the so-called true canopy pixels because this threshold yielded the best results through the visual inspection. This method was based on the approach adopted by Heidarian Dehkordi et al. [26]. Finally, the created WDVI-based mask was applied to the vegetation indices (Table 2) to remove the soil pixels from the analysis. Unsupervised k-means clustering [50] was applied to the multispectral OSAVI imagery in order to generate a MSP crop health map with "good", "moderate", and "poor" classes (Figure 2). The kmeans algorithm identifies the k number of cluster centroids (i.e., 3 in the present study) within the The aforementioned vegetation indices were used for identifying the crop spectral differences between reference and century-old biochar plots throughout the growing season. First, WDVI was used to generate a multispectral-based vegetation mask. A threshold of WDVI higher than 0.30 was used to retrieve the so-called true canopy pixels because this threshold yielded the best results through the visual inspection. This method was based on the approach adopted by Heidarian Dehkordi et al. [26]. Finally, the created WDVI-based mask was applied to the vegetation indices (Table 2) to remove the soil pixels from the analysis.

image and allocates each pixel to the nearest cluster by minimizing its Euclidean distance from the cluster centroid; this iterative algorithm applies until a steady state with the optimized centroid Unsupervised k-means clustering [50] was applied to the multispectral OSAVI imagery in order to generate a MSP crop health map with "good", "moderate", and "poor" classes (Figure 2). The k-means algorithm identifies the k number of cluster centroids (i.e., 3 in the present study) within the image and allocates each pixel to the nearest cluster by minimizing its Euclidean distance from the cluster centroid; this iterative algorithm applies until a steady state with the optimized centroid positions is reached [50]. Moreover, principal component analysis [51] was performed on the red and green spectral channels retrieved from the RGB imagery (Figure 2). The first component (PC1) of the principal component analysis (PCA), which accounts for the maximum proportion of variance [52] from the red and green

bands of the RGB orthomosaic (i.e. 99.06% in the present study), was used in the unsupervised k-means clustering to create a RGB crop health map (Figure 2). More precisely, as recommended by Heidarian Dehkordi et al. [26], the need to explore the impact of century-old biochar on crop health during the senescence phase was addressed in this study.

Moreover, the use of multispectral vegetation indices for the remote sensing-based estimation of crop yield at the field-scale, as suggested by [53], was evaluated in this study (Figure 2). In addition, the impact of century-old biochar on the predicted multispectral crop yield was studied.


**Table 2.** List of the multispectral vegetation indices we have examined in this study.

#### *2.3. Statistical Analysis*

Statistical paired and global t-tests [59] were used to determine whether the contrasts between reference and century-old biochar plots were statistically significant. Several levels of significance were considered to investigate the observed differences. Throughout the manuscript, asterisks \*, \*\*, \*\*\*, and \*\*\*\* indicate the significance levels of *p*-values less than 0.05, 0.01, 0.001, and 0.0001, respectively.

Relationships between UAV-based derivatives and ground-based measurements were evaluated using the statistical F-test. Subsequently, the significance of the obtained correlations was evaluated through the non-parametric Kendall's rank correlation tau (τ) test which is less sensitive to the small statistical sample size, i.e., 11 experimental plots in the present study [60].

#### **3. Results and Discussion**

#### *3.1. UAV-Based Crop Assessment*

#### 3.1.1. Canopy Cover

The temporal pattern of canopy cover shows vegetation development from 20 March to 29 April (Figure 4). At the beginning of May the canopy cover was stable around the maximum until mid-May. From this time onwards, canopy cover slightly decreased towards the end of the season since leaves lost their greenness and began yellowing. The canopy cover clearly showed larger values in biochar plots than the reference plots for most of the season (Figure 4). The Area under the curve (AUC) of the canopy cover was significantly higher for the century-old biochar plots than the reference plots over the season (*p*-value of 0.00007 and 0.00002 for paired and global t-tests, respectively). This finding is in accordance with the positive impact of century-old biochar patches on the evolution of chicory canopy cover in the same experimental farm [26]. Another study by Kerré et al. [6] also identified this positive impact of century-old biochar from relict charcoal hearths on maize crop growth in Belgium.

The greater canopy cover development over the century-old biochar plots could be explained by improved soil physical properties [6] such as water holding capacity (WHC), and soil porosity and infiltration rates, leading to higher water content available to plants that, in turn, increased the vegetation development in the biochar plots. Moreover, higher soil chemical properties (i.e., organic carbon and nitrogen content) were observed over the same century-old biochar plots [26] characterized by a positive relationship with the evolution of canopy cover in the 2018 growing season.

maize crop growth in Belgium.

(AUC) of the canopy cover was significantly higher for the century-old biochar plots than the reference plots over the season (*p*-value of 0.00007 and 0.00002 for paired and global t-tests, respectively). This finding is in accordance with the positive impact of century-old biochar patches on the evolution of chicory canopy cover in the same experimental farm [26]. Another study by Kerré

**Figure 4.** (**a**) Temporal profiles of the canopy cover, derived from the RGB imagery, comparing the 11 century-old biochar plots with the corresponding reference plots throughout the season. (**b**) Box plot visualization of the area under the curve of canopy cover between the century-old biochar and reference plots (i.e., obtained from graphs displayed in sub-plot a for each site). The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. The black circles show all of the experimental plots and the grey dark green line indicates the corresponding pairs of a reference and century-old biochar plot. **Figure 4.** (**a**) Temporal profiles of the canopy cover, derived from the RGB imagery, comparing the 11 century-old biochar plots with the corresponding reference plots throughout the season. (**b**) Box plot visualization of the area under the curve of canopy cover between the century-old biochar and reference plots (i.e., obtained from graphs displayed in sub-plot a for each site). The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. The black circles show all of the experimental plots and the grey dark green line indicates the corresponding pairs of a reference and century-old biochar plot.

The greater canopy cover development over the century-old biochar plots could be explained by improved soil physical properties [6] such as water holding capacity (WHC), and soil porosity and infiltration rates, leading to higher water content available to plants that, in turn, increased the vegetation development in the biochar plots. Moreover, higher soil chemical properties (i.e., organic carbon and nitrogen content) were observed over the same century-old biochar plots [26] characterized by a positive relationship with the evolution of canopy cover in the 2018 growing In addition, the earlier germination of winter wheat plants over the century-old biochar plots on the first acquisition (Figure 4) could be attributed to the dark background soils increasing the soil temperature and subsequently also plant temperature, leading to a faster crop phenological development. Hence, future research using high-resolution thermal imagery should investigate the impact of century-old biochar on soil–plant water exchanges and their relationship with crop development.

#### season. In addition, the earlier germination of winter wheat plants over the century-old biochar plots on 3.1.2. Plant Height

the first acquisition (Figure 4) could be attributed to the dark background soils increasing the soil Figure 5 illustrates the outcome of the height analysis. Winter wheat height increased from 20 March to 29 April (Figure 5) which is in accordance with our earlier finding of canopy cover evolution (Figure 4). The height estimates were stable in May and June around the maximum. Wheat height was rather accurately derived with an accuracy (SD/mean) of 4.92 cm counting the entire dataset (Figure 5), while the previous work by ten Harkel et al. [61] reached a higher accuracy of 3.4 cm for wheat height estimates using UAV-based light detection and ranging (LiDAR) data. The usage of UAV-based LiDAR

development.

3.1.2. Plant Height

enables a proper selection of the top of the canopy pixels to be used in height computations [61] increasing the modeling accuracy, which was not possible in the present study. The highest variation in height estimates (Figure 5) was observed on 28 March (SD of 10.73 and 11.25 cm for the century-old biochar and reference plots, respectively), which could be attributed to the erectophile and open structure of winter wheat at the beginning of the season, lowering the number of true canopy pixels used for height estimates. However, the observed variations were much lower during the maturity phase (with a maximum SD of 5.45 cm on 13 May) because wheat canopies are rather dense, increasing the number of true canopy pixels in height computations. (Figure 5), while the previous work by ten Harkel et al. [61] reached a higher accuracy of 3.4 cm for wheat height estimates using UAV-based light detection and ranging (LiDAR) data. The usage of UAV-based LiDAR enables a proper selection of the top of the canopy pixels to be used in height computations [61] increasing the modeling accuracy, which was not possible in the present study. The highest variation in height estimates (Figure 5) was observed on 28 March (SD of 10.73 and 11.25 cm for the century-old biochar and reference plots, respectively), which could be attributed to the erectophile and open structure of winter wheat at the beginning of the season, lowering the number of true canopy pixels used for height estimates. However, the observed variations were much lower during the maturity phase (with a maximum SD of 5.45 cm on 13 May) because wheat canopies are rather dense, increasing the number of true canopy pixels in height computations.

was rather accurately derived with an accuracy (SD/mean) of 4.92 cm counting the entire dataset

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 9 of 24

temperature and subsequently also plant temperature, leading to a faster crop phenological development. Hence, future research using high-resolution thermal imagery should investigate the impact of century-old biochar on soil–plant water exchanges and their relationship with crop

Figure 5 illustrates the outcome of the height analysis. Winter wheat height increased from 20 March to 29 April (Figure 5) which is in accordance with our earlier finding of canopy cover evolution

**Figure 5.** Comparison of the plant height of all the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median height, respectively. The error-bar presents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum height, respectively. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance; the acronym NS stands for nonsignificant. **Figure 5.** Comparison of the plant height of all the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median height, respectively. The error-bar presents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum height, respectively. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance; the acronym NS stands for non-significant.

The impact of century-old biochar on winter wheat height considering all of the 11 biochar plots versus the reference plots is shown in Figure 5 for each acquisition date. The mean height (mean value of the 11 plots) was higher in the century-old biochar plots than in the reference plots for each acquisition date except 28 March (Figure 5). It is worth noting that the observed mean height differences between reference and century-old biochar plots were statistically non-significant for all the monitoring dates (Figure 5). This slight promotion of wheat height as a consequence of biochar presence was previously reported for lettuce–cabbage–lettuce [62] and oat [63] fields. The greater plant height over the century-old biochar patches could be attributed to the higher soil organic carbon and nitrogen contents as previously found by Schulz et al. [63].

Figure 6 shows the outcome of the height analysis for each experimental plot. Most of the experimental plots showed greater height values in century-old biochar plots compared to their corresponding reference plots throughout the season, except for plots 8 and 11 where an inverse trend was observed (Figure 6). This can be explained by the location of the experimental plots 8 and 11 which were situated in a part of the field with a high soil moisture content (images not shown). The observed

3.1.3. Vegetation Indices

difference of plant height between the century-old biochar plots and the corresponding reference plots at each acquisition date was significant for all of the experimental plots (Figure 6) except for plots 5, 8, and 10 on 28 March and plot 10 on 16 April. The observed height differences between reference and century-old biochar plots were significant for most of the experimental pairs separately (as shown in Figure 6), whereas the differences were not of particular concern when counting the entire dataset (Figure 5). The latter could be explained by the within-field spatial variability of plant height that accounts for the spatial variation in field topographic derivatives such as elevation, groundwater flow, or surface wetness. Validation with the ground-based measurements of plant height should be performed in future studies to explore the robustness of the proposed plant height derivation method in this research. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 11 of 24

**Figure 6.** Box plot visualization of winter wheat height development for the 11 experimental pairs of century-old biochar versus reference plots. The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for nonsignificant. The bottom magnified window displays an example of a height comparison between century-old biochar and reference plot 5 (i.e., on 16 April 2019). **Figure 6.** Box plot visualization of winter wheat height development for the 11 experimental pairs of century-old biochar versus reference plots. The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for non-significant. The bottom magnified window displays an example of a height comparison between century-old biochar and reference plot 5 (i.e., on 16 April 2019).

is completely in line with the trend observed from the canopy cover evolution in Figure 4.

Figure 7 illustrates the temporal profiles for OSAVI. There is a slight vegetation development from 20 March until the end of that month. The OSAVI then increased gradually until the end of

#### 3.1.3. Vegetation Indices

Figure 7 illustrates the temporal profiles for OSAVI. There is a slight vegetation development from 20 March until the end of that month. The OSAVI then increased gradually until the end of April and remained constant at the maximum towards mid-May. From this time onwards, the OSAVI started to decrease as the leaves started to lose their greenness. The overall OSAVI pattern (Figure 7) is completely in line with the trend observed from the canopy cover evolution in Figure 4. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 12 of 24

**Figure 7.** Comparison of the optimized soil adjusted vegetation index (OSAVI) of all of the 11 centuryold biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median OSAVI, respectively. The error-bar represents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum OSAVI, respectively. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for non-significant. **Figure 7.** Comparison of the optimized soil adjusted vegetation index (OSAVI) of all of the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median OSAVI, respectively. The error-bar represents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum OSAVI, respectively. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for non-significant.

The impact of century-old biochar on OSAVI considering all of the 11 century-old biochar plots versus the reference plots for each acquisition date is shown in Figure 7. The mean OSAVI (mean value of the 11 plots) was higher in the century-old biochar plots than the reference plots at each acquisition date except on 24 June when the mean OSAVI was 0.73 for both reference and centuryold biochar plots (Figure 7). This finding is in accordance with the decreased plant greenness at the end of growing season for chicory in century-old biochar patches of the same experimental farm [26]. The observed OSAVI difference was non-significant between the reference and century-old biochar plots towards the end of the season, while it was highly significant from the beginning of the season until mid-April (Figure 7). This could be partially explained by the sparseness of the wheat canopies at the beginning of the season, increasing the interactions between plants and dark background soils in the century-old biochar plots. The impact of century-old biochar on OSAVI considering all of the 11 century-old biochar plots versus the reference plots for each acquisition date is shown in Figure 7. The mean OSAVI (mean value of the 11 plots) was higher in the century-old biochar plots than the reference plots at each acquisition date except on 24 June when the mean OSAVI was 0.73 for both reference and century-old biochar plots (Figure 7). This finding is in accordance with the decreased plant greenness at the end of growing season for chicory in century-old biochar patches of the same experimental farm [26]. The observed OSAVI difference was non-significant between the reference and century-old biochar plots towards the end of the season, while it was highly significant from the beginning of the season until mid-April (Figure 7). This could be partially explained by the sparseness of the wheat canopies at the beginning of the season, increasing the interactions between plants and dark background soils in the century-old biochar plots.

Contrast between reference and century-old biochar plots for each vegetation index is presented in Table 3 for all the acquisition dates. Century-old biochar significantly affected the crop spectral information derived from NDVI, NDRE, OSAVI, CI-red, and s-CCCI on 20 and 28 March. However, Contrast between reference and century-old biochar plots for each vegetation index is presented in Table 3 for all the acquisition dates. Century-old biochar significantly affected the crop spectral

s-CCCI (Table 3). On 29 April only NDVI, NDRE, and CI-red, and on 13 May only NDVI and CI-red, exhibited a significant impact of century-old biochar on crop spectral status. All of the examined vegetation indices exhibited non-significant differences between reference and century-old biochar plots on 24 June (Table 3). This finding is in accordance with [26] who reported no significant information derived from NDVI, NDRE, OSAVI, CI-red, and s-CCCI on 20 and 28 March. However, WDVI, CVI, and EVI did not exhibit significant differences in crop spectral status between reference and century-old biochar plots during March. All of the examined vegetation indices showed a significant impact of century-old biochar on crop spectral status on 16 April except WDVI, EVI, and s-CCCI (Table 3). On 29 April only NDVI, NDRE, and CI-red, and on 13 May only NDVI and CI-red, exhibited a significant impact of century-old biochar on crop spectral status. All of the examined vegetation indices exhibited non-significant differences between reference and century-old biochar plots on 24 June (Table 3). This finding is in accordance with [26] who reported no significant difference of the multispectral vegetation indices between the century-old biochar and reference plots for chicory crop at the beginning of July.

**Table 3.** Mean vegetation index values over the study area comparing the century-old biochar and reference plots during the growing season of 2019. The statistical significance of the difference is expressed as the p-value of the paired t-test. Next to the p-values asterisks, \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and acronym NS stands for non-significant.


It is worth mentioning that the WDVI and the EVI performed the worst in discriminating spectral difference between century-old biochar and reference plots throughout the season. The latter could be attributed to the use of the blue band in EVI's computation which is not of particular importance in the vegetation spectral signature [64].

Figure 8 shows the outcome of the OSAVI analysis for each experimental plot. Most of the experimental plots exhibited higher OSAVI values in century-old biochar plots compared to their adjacent reference plots for most of the season (Figure 8). Graphs of the other examined multispectral vegetation indices (Table 2) are not presented here. No clear trend was observed in terms of the OSAVI difference between reference and century-old biochar plots (Figure 8). The contrast between reference and century-old biochar plots was highly significant at several acquisitions (e.g., for plot 1 on 29 April), whereas some acquisitions exhibited a non-significant difference (e.g., plot 1 on 28 March). *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 14 of 24

**Figure 8.** Box plot visualization of optimized soil adjusted vegetation index for the 11 experimental pairs of century-old biochar versus reference plots over the growing season of 2019. The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for non-significant. No spectral information is available for reference plot 1 on 20 March because of a flight planning constraint. **Figure 8.** Box plot visualization of optimized soil adjusted vegetation index for the 11 experimental pairs of century-old biochar versus reference plots over the growing season of 2019. The horizontal black line displays the median value, surrounded by box edges representing the 25th and 75th percentiles. Outside of the plot, asterisks \*, \*\*, \*\*\*, and \*\*\*\* reveal the statistical levels of significance and the acronym NS stands for non-significant. No spectral information is available for reference plot 1 on 20 March because of a flight planning constraint.

#### 3.1.4. Crop Health 3.1.4. Crop Health

retrieve an accurate crop health map.

Both canopy cover and OSAVI were characterized by a decreasing trend on 24 June (Figures 4 and 7, respectively). The impact of century-old biochar on crop health during the senescence phase is presented in Figure 9. In addition, Figure 9 presents the two different crop health maps, obtained from the two different sensors, including the associated general methodological workflows as represented in Figure 9a,b for the multispectral sensor and Figure 9c,d for the RGB sensor. An unsupervised k-means clustered crop health map of OSAVI including "good", "moderate", and "poor" classes, corresponding to 24 June, is presented in Figure 9b. In addition, an unsupervised kmeans clustered crop health map derived from the first component (PC1) from the principal component analysis of the red and green bands of the RGB orthomosaic image corresponding to 24 June is shown in Figure 9d, clustering the study area into "good", "moderate", and "poor" classes. It is worth mentioning that 99.06% of the observed variance was explained by the first component (PC1) and only 0.94% of the variance retained by the second component (PC2) (Figure 9c). Moreover, there Both canopy cover and OSAVI were characterized by a decreasing trend on 24 June (Figures 4 and 7, respectively). The impact of century-old biochar on crop health during the senescence phase is presented in Figure 9. In addition, Figure 9 presents the two different crop health maps, obtained from the two different sensors, including the associated general methodological workflows as represented in Figure 9a,b for the multispectral sensor and Figure 9c,d for the RGB sensor. An unsupervised k-means clustered crop health map of OSAVI including "good", "moderate", and "poor" classes, corresponding to 24 June, is presented in Figure 9b. In addition, an unsupervised k-means clustered crop health map derived from the first component (PC1) from the principal component analysis of the red and green bands of the RGB orthomosaic image corresponding to 24 June is shown in Figure 9d, clustering the study area into "good", "moderate", and "poor" classes. It is worth mentioning that 99.06% of the observed variance was explained by the first component (PC1) and only 0.94% of the variance retained by the second component (PC2) (Figure 9c). Moreover, there was a positive significant correlation between the red and green bands (Kendall's τ coefficient = 0.89 and *p*-value(Kendall) < 0.0001). The experimental field was completely covered by vegetation on 24 June enabling a proper PCA analysis over the vegetation pixels, without any soil background effects, to retrieve an accurate crop health map. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 15 of 24

**Figure 9.** Data from 24 June 2019: (**a**) Map of optimized soil adjusted vegetation index (OSAVI). (**b**) Crop health map derived from the k-means clustering of OSAVI; the magnified window represents an example of an area characterized by poor crop health in the study field. (**c**) RGB orthomosaic image (left), the green and red spectral responses of the visible RGB orthomosaic image (middle), and the first component (PC1) raster of the principal component analysis (PCA) of the green and red bands (right). (**d**) Crop health map derived based on the k-means clustering of the PC1 raster of the PCA of green and red spectral channels. **Figure 9.** Data from 24 June 2019: (**a**) Map of optimized soil adjusted vegetation index (OSAVI). (**b**) Crop health map derived from the k-means clustering of OSAVI; the magnified window represents an example of an area characterized by poor crop health in the study field. (**c**) RGB orthomosaic image (left), the green and red spectral responses of the visible RGB orthomosaic image (middle), and the first component (PC1) raster of the principal component analysis (PCA) of the green and red bands (right). (**d**) Crop health map derived based on the k-means clustering of the PC1 raster of the PCA of green and red spectral channels.

Categorical classification of the experimental plots into the crop health clusters is presented in Table 4 for both of the multispectral and RGB crop health maps (Figure 9). Based on the multispectral crop health map, the average health was 42.5% for the century-old biochar plots compared to 39.4% in the reference plots. Similarly, the average health based on the RGB crop health map was higher for the century-old biochar plots (54.6%) than the reference plots (51.6%). This finding shows that there is only a small positive impact of century-old biochar on crop health status during the senescence phase. The present study paves the way for further research focusing on the effects of biochar on crop health during the senescence phase. Categorical classification of the experimental plots into the crop health clusters is presented in Table 4 for both of the multispectral and RGB crop health maps (Figure 9). Based on the multispectral crop health map, the average health was 42.5% for the century-old biochar plots compared to 39.4% in the reference plots. Similarly, the average health based on the RGB crop health map was higher for the century-old biochar plots (54.6%) than the reference plots (51.6%). This finding shows that there is only a small positive impact of century-old biochar on crop health status during the senescence phase. The present study paves the way for further research focusing on the effects of biochar on crop health during the senescence phase.


**Table 4.** Categorical classification of the crop health into good, moderate, and poor classes comparing the century-old biochar versus reference plots. The classification is based on the multispectral (MSP) crop health (left) and the first component (PC1) of the principal component analysis (PCA) of the green and red spectral channels of the RGB orthomosaic image (right).

The cross-tabulation analysis of the computed crop health maps is presented in Table 5 comparing the number of all of the pixels within the study area classified into good, moderate, and poor health classes. The result yielded a clustering agreement of 74.82% (Table 5). Out of a total of 564,760 pixels across the entire study area, 36.06% were classified as good, 31.54% as moderate, and 7.21% as poor by both methods (Table 5). Ground-based optical imagery, such as from the Fabry-Perot interferometer (FPI) camera system [39], is an essential step for future research focusing on the evaluation of crop health because of its improved spatial and spectral resolutions, enabling an accurate monitoring of the development of various plant diseases at the leaf level.

**Table 5.** Cross-tabulation table of the computed MSP and RGB crop health maps comparing the clustering accuracy of all of the pixels within the study area into good, moderate, and poor classes.

