2.2.2. UAV-LiDAR Data

LiDAR data were collected using the RIEGL RiCOPTER UAV system (Riegl, Austria) with the VUX-SYS laser scanner mounted underneath. The authors of [15] provided a description of the details of the system. The VUX LiDAR unit has an accuracy of 1 cm and precision of 0.5 cm. The absolute positioning error on the ground also depends on the IMU accuracy, which was <5 cm in the horizontal and <10 cm in vertical direction. Analysis of previous datasets acquired with the same setup resulted in positioning errors of <5 cm. The scanner pulse repetition rate was set at 550 kHz, the scanner angle from 30–330 degrees, and the scanner speeds were synchronized with the UAV forward speed to create a regular point spacing. On the predefined dates, UAV-LiDAR data was acquired by flying at 40 m above ground level (m.a.g.l.) with a programmed speed of 6 m/s. The actual flight altitude and speed can deviate slightly from those preset values and were calculated from the flight recordings afterwards (Table 4: mean flying height and speed).

Further, to assess the influence of the flying height and speed on the accuracy of acquired phenological parameters, additional flights were done within a few days after the second flight day, with heights varying from 20 to 90 m above ground level, and with programmed speeds ranging from

3 m/s to 8 m/s. An overview of the flights and the actual flight data is shown in Table 4. The code given in Table 4 is later used to identify the data acquired from each flight.

**Table 4.** Overview of flight dates, including flight codes and corresponding general flight details The flight codes are 1) indicating the flight number; 2) the flight altitude (LA = Low Altitude, MA = Medium Altitude, HA = High Altitude; 3) and flight speed (LS = Low Speed, MS = Medium Speed, HS = High Speed). Those codes will be used throughout the paper.


If during the analysis it appeared that field measurements of plant height or biomass ended up in the same pixel as a driving path as seen from the Lidar 3D point cloud, the point was moved 1 pixel away from the driving path. This was the case for one point in the winter wheat dataset. Further, the field measurement locations were selected before the flights were made. Points that appeared to be in areas with a much lower LiDAR point density were excluded from the height and biomass analysis. Two winter wheat points during the last flight date had to be excluded for this reason.

#### *2.3. Data Pre-Processing*

Riegl RiCOPTER data was pre-processed to co-registered point cloud datasets (Figure 2) using the standard Riegl processing chain as outlined by Brede et al. [15]. This includes trajectory processing, for which a virtual GNSS base station was used. The movements of the UAV platform were corrected using the POSPac software (Applanix, Canada), which combines the GNSS correction and the movement data recorded by the IMU, with a reported precision of ~1–2 cm for all flights. To produce a geo-referenced point cloud, the flight path, and raw point cloud data were combined in RiPROCESS (Riegl, Austria). Further refinement of the point cloud, using automatically detected surfaces, was done by RiPRECISION (Riegl, Austria). The final co-registered and geo-referenced point cloud datasets were stored in LAS format.

**Figure 2.** Overview of methodology for height and biomass retrieval from UAV-LiDAR observations for crops.

The point clouds were clipped for the boundaries of the three fields (Figure 1). The ground points were selected using LasTools (Rapidlasso, USA) with the lasground function. From the points classified as ground, a triangulated irregular network (TIN) was created. The height of each non-ground classified point was calculated as the height above this TIN, using the "replace z" option in lasground. The same procedure was applied for all the three fields.

#### *2.4. Data Analysis*

#### 2.4.1. LiDAR Based Plant Height

Due to the difference in the structure of the crop canopies, the highest point in the point cloud is not always the best representation for the field-measured crop height. Therefore, for each grid cell of 1 m2 the average plant height was derived by including a certain number of top-of canopy points. To determine the amount of top-of-canopy points needed to get the best plant height estimation, a theory was used that is more commonly used in economics to decide on the best allocation of products: the Pareto efficiency score [16]. This method compares different alternatives and based on preset criteria, calculates a score of optimal performance. In this case, we used this method here to decide based on two criteria, namely the highest R<sup>2</sup> and lowest root mean squared error (RMSE), which are calculated for each alternative of measured plant height and a certain number of top-of canopy points. The method then only selects the combination for which no better combination exists. This can result in multiple optimal Pareto efficiency scores. Next, from these Pareto efficient combinations, the combination with the highest R<sup>2</sup> is selected. This was done once for the whole growing season for every crop. The corresponding number of top-of canopy points of this combination was used to calculate the plant height of each raster cell within the field. A linear model relating field plant height measurements to LiDAR-derived plant height was built using the 3D LiDAR point cloud by taking the optimal number of top-of-canopy points into account.

#### 2.4.2. Crop Biomass

To estimate biomass from the LiDAR point cloud, the 3DPI method as developed by the authors of [2] and described by Equation (1) was used and calculated for a cell size of 1 m2:

$$\text{GDPI} = \sum\_{i=0}^{i=n} \left( \frac{p\_i}{p\_t} e^{k \frac{p\_{ci}}{p\_t}} \right) \tag{1}$$

where *i* is a given 10 mm vertical layer with 0 being the ground layer and n the uppermost layer respectively; *k* is a tuning parameter changing from −1 to 5 by steps of 0.05; *pi* is the number of LiDAR points for a given layer of 50 mm; *pt* is the total number of LiDAR points for all layers; and *pcs* is the cumulative sum of LiDAR points intercepted above a given layer of 50 mm, as described in [2].

To derive the biomass prediction model, a linear regression was performed between the dimensionless 3DPI indicator and the biomass measured in the field based on the 3DPI indicator. First, the optimal value for the parameter k was determined. For each individual crop and flight, the 3DPI score was calculated by varying the k-parameter from −1 to 5 by increments of 0.05. The best k-parameter was chosen based on the correlation between the field measurements and the 3DPI score. The data for each crop was combined into a crop-specific seasonal prediction model for the flights Day1-MA-MS, Day2-MA-MS, and DAY4-MA-MS. Secondly, all the data of these flights and crops were combined into a general non-crop-specific seasonal prediction model.

Next, the flight data from DAY2-LA-MS, DAY2-LA-LS, DAY3-HA-HS, and DAY3-LA-LS were used to determine the effects of different flight specification and were compared to Day2-MA-MS. The 3DPI for these four flights was calculated using the same field measurements of plant height and biomass as for Day2-MA-MS. This enabled the analysis of different flight settings because only the 3DPI indicator differed between the flights so a direct comparison can be made. To visualize the

differences between these four flights, maps of predicted biomass were derived and difference maps of each flight with flight Day2-MA-MS were calculated.

#### 2.4.3. Model Validation

To assess the performance of the models (Figure 2) for determining biomass, plant height, and the general models, three statistical indicators were used; the coefficient of determination (R2), the mean absolute error (MAE), and the RMSE. We included a cross validated MAE and RMSE, using a k-fold cross validation. For each k-fold, the out-sample RMSE and MAE were calculated, after which they were averaged over the n-number of k-folds used in the cross validation. These averaged out-sample RMSE and MAE were compared to the in-sample RMSE and MAE. This comparison shows the predictive power of the model (out-sample errors).

For the plant height analysis, a 10-fold cross validation was performed, using 61, 51, and 51 samples for potato, sugar beet, and wheat, respectively (Table 3). For the biomass estimation, a repeated 5 × 3-fold cross validation was performed. The reason for this repeated cross validation for biomass is the limited amount of biomass samples available for the statistical analysis. In this case, we used 15, 15, and 9 samples for potato, sugar beet, and wheat, respectively (Table 3).

For the biomass error analysis, the MAE and RMSE were then normalized by dividing by Xmax-Xmin, where Xmax is the highest measured biomass sample in the dataset and Xmin the lowest measured biomass sample in the dataset, resulting in a normalized mean absolute error (NMAE) and normalized RMSE (NRMSE).
