Deriving Aerodynamic Roughness Length at Ultra-High Resolution in Agricultural Areas Using UAV-Borne LiDAR

Abstract

The aerodynamic roughness length (Z₀) and surface geometry at ultra-high resolution in precision agriculture and agroforestry have substantial potential to improve aerodynamic process modeling for sustainable farming practices and recreational activities. We explored the potential of unmanned aerial vehicle (UAV)-borne LiDAR systems to provide Z₀ maps with the level of spatiotemporal resolution demanded by precision agriculture by generating the 3D structure of vegetated surfaces and linking the derived geometry with morphometric roughness models. We evaluated the performance of three filtering algorithms to segment the LiDAR-derived point clouds into vegetation and ground points in order to obtain the vegetation height metrics and density at a 0.10 m resolution. The effectiveness of three morphometric models to determine the Z₀ maps of Danish cropland and the surrounding evergreen trees was assessed by comparing the results with corresponding Z₀ values from a nearby eddy covariance tower (Z₀_EC). A morphological filter performed satisfactorily over a homogeneous surface, whereas the progressive triangulated irregular network densification algorithm produced fewer errors with a heterogeneous surface. Z₀ from UAV-LiDAR-driven models converged with Z₀_EC at the source area scale. The Raupach roughness model appropriately simulated temporal variations in Z₀ conditioned by vertical and horizontal vegetation density. The Z₀ calculated as a fraction of vegetation height or as a function of vegetation height variability resulted in greater differences with the Z₀_EC. Deriving Z₀ in this manner could be highly useful in the context of surface energy balance and wind profile estimations for micrometeorological, hydrologic, and ecologic applications in similar sites.

1. Introduction

To achieve more sustainable land and water-use management in precision agriculture, the serviceable description of the complex interactions between the vegetation growth, water availability, and aerodynamic characteristics of a canopy is pivotal. In this context, aerodynamic roughness length for momentum transfer (Z₀) is a key factor in micrometeorological and hydrological applications, since it can elucidate how surface geometry may lead to alterations in energy, gas, and water exchanges surface friction; and the deflection of airflow [1,2,3]. Many models have been developed to estimate the components of surface energy balance and evapotranspiration (ET) [4,5,6,7] using passive remote sensing observations and a set of algorithms to retrieve surface parameters such as the aerodynamic resistance to heat transfer (r_ah), which is a function of Z₀ [8,9]. Obtaining an efficient parameterization for r_ah has been a challenging task, and there is no single method to accurately estimate Z₀ over a wide range of land cover types [10,11], introducing further uncertainty in the modeling of energy fluxes and ET. As such, incremental advances in resolving Z₀ using canopy structure data at ultra-high resolution from drone-borne instrumentation have the potential to improve the accuracy of surface energy balance models in the context of energy, gas, and water exchange estimations in precision agriculture and related applications.

Two primary approaches exist in assigning Z₀ values. The most commonly used method is based on micrometeorological observations obtained by an eddy covariance (EC) system and the Monin–Obukhov similarity theory. A disadvantage of this approach is that Z₀ is restricted to single average values in the flux footprint of the EC system. Moreover, estimates of Z₀ cannot always rely on field-based experiments due to practical limitations and high costs. The second approach relies on the demarcation of surface objects using remote sensing observations and the establishment of empirical relationships between Z₀ and measurable characteristics of site-specific roughness elements. In this morphometric method, theoretical models of the boundary layer are combined with more sophisticated physical models of vegetation canopy to determine Z₀ (e.g., [12]). Following this approach, Z₀ is often associated with the frontal area index (fai), which is the area of windward vertical faces of the roughness elements to the total area under consideration, and the plan area index (pai), which equals the horizontal area occupied by roughness elements divided by the total area [13]. For agricultural and natural sites without density information, Z₀ is often simply related to canopy height (h) [14]. However, this relation is not always constant, since the density, the type of vegetation, and the micro- or macrotopographic characteristics can affect Z₀ variations as well [15].

Using remote sensing optical imagery and morphometric models, numerous researchers have developed semi-empirical formulas applicable to obtain fractional values of Z₀/h of vegetation surfaces [16], and some of them have included the effects of vegetation indices (VIs) [17]. However, precise observations of vegetation height or VI in both high spatial and temporal resolution are difficult to obtain from publically available satellite datasets [18]. A major weakness of airborne/satellite imagery is its limitation in viewing beneath the canopy, leading to sparse points and low-density information on bare soil [19], whereas light detection and ranging (LiDAR) scanners can provide quantitative information on the 3D structure of a canopy because the laser pulses can partly penetrate vegetation cover. The technology of airborne LiDAR scanners (ALS) has now successfully been employed for the extraction of surface roughness characteristics in forests [20,21], urban areas [22], and low-vegetation areas [23]. The representation of a mixed grassland prairie by ALS datasets revealed that up to 76% of the variation in Z₀ was due to the height variability of vegetation and up to 65% of the variation could be explained by estimates of vegetation height [24]. Li et al. [25] found that the accuracy of the estimated Z₀ of a semi-arid shrubland using ALS data depends on the adopted morphometric models and the precise representation of shrub height in these models. In short and dense canopies, the estimation of vegetation height using ALS is prone to errors [26], mainly due to the lack of identifiable referenced objects and of detectable differences between first (i.e., vegetation) and last (i.e., ground) LiDAR returns [27,28].

Given these challenges, the advantages conferred by unmanned aerial vehicle (UAV)-borne-LiDAR scanners may yield fine-grained spatiotemporal estimations of Z₀ in agricultural areas. Compared to manned ALS, the comparatively cost-effective UAV-LiDAR systems are more flexible in data sampling and produce higher point cloud density due to the larger field of view of the scanner and lower flight altitude and speed, allowing a larger number of LiDAR beams per scan [29]. These characteristics may limit the commonly observed underestimation of canopy heights and mitigate difficulties in deriving individual roughness element canopies from airborne LiDAR data [30]. Resop et al. [31] documented that higher-resolution UAV-LiDAR data facilitated the identification of small vegetation and micro-alterations in a heterogeneous terrain that were not detectable by ALS observations. In a similar study, the 3D characterization of individual plant species of a shrubland area was achievable at the submeter scale using a UAV-LiDAR system [32]. However, the technology of UAV-LiDAR is not currently used in precision farming, despite its ability to effectively monitor canopy density [33] and fine-scale variations in crops attributes compared to UAV-optical imagery [34,35,36], which is widely employed in such applications [37,38,39,40].

Most published methods for the estimation of crop height from LiDAR sensors are based on the determination of the canopy height model (CHM), which is derived from the difference between rasterized ground elevation data (digital terrain model (DTM) and rasterized original point cloud elevation data (digital surface model (DSM)). A DTM of an agricultural field can be easily retrieved by scanning the bare soil before vegetation growth, but such data acquisitions are not always feasible due to land management practices, so a variety of filtering methods to generate DTMs and CHMs from airborne LiDAR data have been developed. As automatic separation of ground and non-ground points from LiDAR data has proved to be challenging and the filtering algorithms typically have problems distinguishing ground returns and points reflecting vegetation [41], the adaptability of existing approaches for processing point cloud obtained by UAV-LiDAR is currently inadequate.

Poor filtering performance, particularly in areas of dense vegetation that hides underlying terrain features, may result in erroneous surface morphologies or in sparse ground points, which, when interpolated, fail to reproduce surface morphology. In this framework, a detailed understanding of the errors associated with producing CHMs using UAV-LiDAR technology and a context-specific approach to assess filtering algorithms and morphometric methods to estimate surface characteristics are required.

We explored the potential of UAV-borne LiDAR in the estimation and mapping of surface roughness of a typical dense agricultural environment at very high spatial resolution. The major components of the project were as follows:

• An evaluation of the performance of three segmentation approaches (i.e., a morphological filter (MF), a progressive triangulated irregular network densification filter (TIN), and a combination of MF and TIN) to reliably partition the UAV-LiDAR-derived point cloud data into bare earth and vegetation and, consequently, to generate CHMs at centimeter resolution.
• An assessment for calculating Z₀ values from UAV-LiDAR data using three different morphometric methods (Kustas et al. [14], Raupach, [42], and Menenti and Ritchie [43]), and a comparison of the respective results with the Z₀ obtained from EC observations.
• A discussion of the challenges and further potential of UAV-LiDAR in precision agriculture and related applications.

2. Materials and Methods

2.1. Site Description

The experimental site is part of the Integrated Carbon Observation System network (ICOS) located in Jutland, Denmark (56.037644° N, 9.159383° E). The site is a dense agricultural area with negligible topographic relief covered by potato plants with heights varying from 0.3 to 1.7 m and sparse trees (Figure 1). The climate is humid temperate; during the summer of 2019, the mean temperature was 18.54 °C and the mean relative humidity was 66.65%. The site’s prevailing winds originated from the west (199°) and the mean wind speed recorded was 2.83 m/s.

2.2. Data Collection

The aerial campaign was conducted on a monthly scale before and during the germination of vegetation (14 May, 26 June, 14 July, and 12 August 2019). Point cloud data were acquired using a UAV-LiDAR system (LidarSwiss GmbH, CH) onboard an octocopter Matrice 600 Pro drone (Figure 2). The LiDAR system included an internal navigation system (INS—Oxts Xnav 550 OEM IMU/GPS system) that fused data from an inertial measurement unit (Oxts micro electro-mechanical systems) and GPS data received by a Global Navigation Satellite System antenna, a beam LiDAR scanner (Quanergy M8), a 20 mp SONY RGB camera (16 mm lens), and an integrated data storage unit. A Trimble real-time kinematic GNSS base station was used to provide additional overhead communication with the INS. The horizontal field of view (FOV) of the laser scanner was 360 degrees and the vertical FOV was 20 degrees. The LiDAR data were recorded at 40 m above the ground, producing an average point density of 250 points/m² (max 400 points/m²). The swath width of a single pass was 89 m, and the overlap between two adjacent swaths was greater than 25%. The raw GNSS data files obtained by the INS were converted to position data in pos format using trajectory software (RT Post-process of the NAVsuite software package). The laser scanner’s data were initially produced in bin format and were converted to point clouds in las format using the position data (Geo-LAS software).

An EC system was installed at the western edge of a fenced measuring plot (80 × 20 m), approximately in the center of the agricultural field (orange point in Figure 1). The instrumentation included a Gill R3-50 sonic anemometer (Gill Instruments Ltd., Lymingdon, UK) and a LI7000 closed-path infrared CO₂/H2O gas analyzer (LI-COR, Lincoln, NE, USA). During the sampling period, EC data were recorded at a nominal sampling frequency of 20 Hz and ancillary meteorological data at 1 Hz (for further details, see [44]).

2.3. Evaluation Procedure

A relative homogeneous subscene area of the agricultural field (100 × 150 m covered area) was used to evaluate the performance of the different filtering algorithms applied to the LiDAR dataset retrieved in June (Plot 1 in Figure 1). The optimal set of parameters for each filter was estimated by comparing their effectiveness in terms of the total errors (TEs) of misclassification. Each filter with the optimized parametrization was then applied to the point cloud data, representing the subscene Plot 1 for July and August, and to a second subscene (100 × 270 m) covered by low and high vegetation (Plot 2 in Figure 1). The manual classification of point clouds into terrain and vegetation cannot entirely be employed in a GIS framework because filtering errors in such dense areas are not so obvious to interpret with the naked eye. Instead, bare earth elevation data, as resourced in May before vegetation growth, were combined with the rest of the collected point cloud data from June to August in order to manually label the non-ground points as referenced data of vegetation. The plant heights, represented by the CHMs produced at 0.10 m resolution, were compared with the plant heights that were manually measured in randomly selected 1 × 1 m geolocated plots during the aerial campaign. Based on the CHMs’ various geometric parameters, we evaluated the effectiveness of the morphometric models to obtain Z₀ values by comparing them with the anemometric-based method for specific flux footprint areas.

2.4. Point Cloud Processing

Due to multi-path reflection, the pulses emitted from a LiDAR scanner may reach the ground without returning directly to the instrument but rather reaching neighboring surface points, creating noise data within the point clouds. By calculating the standard deviation of each point’s surrounding fitting plane and by defining an expected relative error (RE = 2), a point was labeled as an outlier if its distance from the fitting plane was greater than the mean average distance plus the product of relative error and standard deviation [45]. Additionally, low noise points that were close to the ground were excluded from further analysis by comparing the maximum height difference between each point and their neighboring points with a height threshold (max height = 0.2 m).

The classification of the point cloud data to vegetation and bare earth was based on the triangulated irregular network densification filter (TIN) [46], a morphological filter (MF) [47], and progressive triangulated irregular network densification (PTD), introduced by Zhao et al. [48]. The efficiency of a filtering algorithm depends on the choice of opening window/grid step and threshold used within each filter. These values are uncertain in all filters and are expected to depend on the size of the objects and the land cover type [49].

The morphological filter (MF) is based on a progressively increasing window size in an iterative process. If the elevation difference between the original data and the data after the opening operation is higher than a user-defined elevation threshold (zthresh), the grid is labeled as a non-ground grid. A number of runs was conducted for the MF filter using a maximum window size from 0.4 to 1 m with an increment of 0.2 m and an elevation threshold from 0.25 to 1.5 m in 0.25 m intervals using a commercial software package [50].

To identify likely ground points in the TIN method, some of the terrain-local minimum points were used as the initial ground seed points to build an initial triangulated irregular network. The sensitivity of the filter was tested by setting it in a software package [51]: a grid step of 0.2 to 1.6 m in 0.2 m intervals based on practical experience, and the spike parameter to 0.10 m, 0.25 and 0.5 m considering the standard deviation for planar patches equal to 1. The spike parameter describes the distance above the coarsest triangulated network for which the points are classified as terrain.

In the PTD algorithm, ground seed points are acquired through a morphological opening operation instead of using the lowest points in user-defined grids. The parameters of iterative distance were set to vary from 0.2 to 1 m in an increment of 0.2 m, and an iterative angle from 2° to 18° in an increment of 4°, using a commercial software package [52].

To assess the performance of the three applied filtering methods, total error (TE) was calculated according to the following equation [53]:

TE = (a + b)/e

(1)

where a is the number of ground points that have been incorrectly classified as vegetation points, b is the number of vegetation points that have been incorrectly classified as ground points by comparing the processed point cloud datasets to the referenced data for Plot 1, and e is the total number of points tested.

Once vegetation was segmented from the LiDAR dataset, the resultant ground points were interpolated to replace non-ground points with an approximation of the correct surface morphology. The inverse distance weighted method [54] was used to interpolate the voids (Figure 3). The point clouds were converted to raster gridded elevation layers of 0.03–0.05 m resolution by connecting all the available point features into a network of triangles and interpolating over the triangular faces using the feature elevation and slope values. To reduce the noise in the raster images by using all the point features, the point cloud data were spatially binned into areas corresponding with the size of the output grid cells [50]. In our analysis, one elevation value from each of the spatial bins was used to generate a gridded layer with a 0.10 m resolution.

2.5. Description of Morphometric Methods

Three morphometric models were assessed to determine the Z₀ of the subscenes, Plot 1 and Plot 2, surrounding the EC tower through surface morphology.

2.5.1. Roughness Length Based on Vegetation Height

The rule of thumb (RT) method only requires the average roughness element height ($\overline{h}$) per pixel, which is linearly related to Z₀ [14] as follows:

Z₀_RT = 0.1 h

(2)

2.5.2. Roughness Length Based on Vegetation Geometry and Wind Conditions

In more advanced morphometric models, the alterations and resulting effects of canopy drag can be included by calculating the drag coefficient [55]. Raupach’s model (RAP), described in Equation (3a, 3b), includes the drag coefficient of an isolated roughness element (c_s = 0.003); the drag coefficient for the substrate surface at $\overline{h}$ (c_r = 0.3); the roughness sublayer influence function (ψ_h = 0.193, accounting for the correction to the logarithmic wind profile); the wind speed (U), the friction velocity, u_*, (u_*/U)max = 0.3; and a free parameter (c_d1 = 7.5). Z_d (m) is the zero-plane displacement and k is the von Karman’s constant (= 0.4).

Z₀_RAP = h (1 − Z_d/h) exp(−k U/u_* + ψ_h)

(3a)

Z_d/h = 1 + {(exp [−(2c_d1 fai)^0.5] − 1)/(2c_d1 fai)^0.5}
u_*/U = min [(c_s + c_r fai)^0.5, (u_*/U)max]

(3b)

The frontal area index (fai) can be defined by integrating positive height changes (Δy) over a cross-sectional line divided by the distance (Δx) in that section length, assuming an isotropic surface [25,56].

fai = (∑Δy)/(∑Δx) for Δy > 0

(4)

The advantage of this technique is that we do not consider the precise shape of a plant but rather its cross-sectional area perpendicular to the wind [22]. In this analysis, we extracted cross-sections along 360/24 sectors reflecting different wind directions from each generated CHM to derive the frontal area indexes. The associated morphometric Z₀ from all directions were averaged into one value at each grid cell.

The plan area index (pai), which equals the horizontal area occupied by roughness elements divided by the total area under consideration, was also calculated since it can be associated with Z₀ and Z_d [15]. For instance, it was observed that as surface cover increases, the magnitude of Z_d/$\overline{h}$ produces a convex curve asymptotically increasing from zero to unity, which is the maximum possible value of pai. The pai and fai for each wind direction were calculated using the UMEP plugin [57] in the open-source geographical information software QGIS [58].

2.5.3. Roughness Length Based on Vegetation Height Variability

This empirical model of Menethi and Ritchie [43] (MR) determines Z₀ as a function of vegetation height variability for each grid cell that is segmented by subcells following:

Z₀_MR = (1/N) ∑ (σi,j/hi,j) h_avg

(5)

where N is the number of subcells within each grid cell; σi,j is the standard deviation of the LiDAR-derived vegetation height (hi,j) per each subcell I; j. h_avg is the average vegetation height calculated from the LiDAR’s CHM. It was documented that coarser grid cells reduce the standard deviation of height regardless of the size of the subcells, while larger subcells lead to higher values of Z₀ [25]. Based on these observations, the size of each grid cell was chosen to be equal to 1 m and the segment size inside each grid was 0.25 m, reflecting the maximum expected variance in plant height within a 1 × 1 m cell. All the geometric parameters of vegetation from the CHMs were retrieved using QGIS.

2.6. Description of the Anemometric Method

The raw data of wind, carbon dioxide, water vapor, and sonic temperature data were processed using EddyPro 4.2.1 software (LI-COR, Lincoln, NE, USA) to estimate half-hourly turbulent scalar fluxes. The processing included statistical tests for raw data screening [59], angle of attack correction [60], 2D coordinate rotation, block averaging, time lag optimization to maximize covariance between vertical wind speed and gas concentrations, humidity corrections applied to the sonic temperature [61], and compensation for density fluctuations [62]. Flux quality flags were assigned according to the test for steady-state conditions and fully developed turbulence following Foken et al. [63] and simplified by CarboEurope-IP.

From the EC micrometeorological data, the estimations of friction velocity, sensible heat flux, and Obukhov length, as well as the meteorological data of wind speed and direction, were used to calculate Z₀ (Z₀_EC) following the logarithmic wind law. For stable or unstable atmospheric conditions, the logarithmic wind profile [64] is given by:

Z₀_EC = (z − Z_d)/exp(kU/u_* + ψ_m)

(6)

where z is the measurement height (m). The stability correction for momentum ${\psi }_m$ was computed using the parameterizations suggested by Högström [65]. Under neutral atmospheric conditions, ${\psi }_m$ equals zero. Z_d was considered here to be equal to 0.7 h [14].

Half-hourly estimations of Z₀_EC from 25 to 27 June, from 13 to 15 July, from 12 to 13 of August, and for daytime hours (from 9:00 to 19:00 local time) were used as reference values for validating the morphometric-derived Z₀. This EC dataset was selected for the comparison analysis of the different methods to calculate Z₀ in order to minimize the effect of the differences in the temporal and spatial resolution of the remote sensing data acquired on the 26 June, 14 July, and 12 August, and in situ EC data.

3. Results

3.1. Segmentation of Point Cloud Data

The morphological filter performed poorly with small window size values when the optimal size of the opening operation was close to 1 m due to the lack of a sufficient number of ground points (Figure 4a). There were fewer Type I errors (accepting a ground point as a vegetation point) than Type II errors (accepting a non-ground point as ground) because there were many local maxima for the filter to identify. When zthresh was increased to 0.5, the Type I errors reduced, meaning that the filter could distinguish morphological differences between ground and non-ground. The optimum threshold was 1 m (

Table 1. Optimized parameters of each filtering algorithm for the tested point cloud dataset.

Methods	Parameter Set
	Window Size (m)	Elevation Threshold (m)	Iterative Distance (m)	Iterative Angle (°)	Grid Size (m)	Spike (m)
MF	1	1
PTD			0.4	4
TIN					0.8	0.2

) and Type II errors increased as zthresh increased beyond 1.5 m since fewer non-ground points were classified correctly. Generally, the filter was not as sensitive to the value of zthresh, performing nearly as well as the optimal parameter combination across a range of values of elevation threshold (Figure 4a). In PTD, the iterative distance is related to the topographic relief that would be expected to be relatively small for flat terrains. From the sensitivity analysis, the iterative distance became stable when it was greater than 0.4 m for an iterative angle equal to 4° (Figure 4b). The TIN model was more sensitive to the grid size, which should be as large as the size of the biggest object located in the filtered area. The suggested step size is between 0.6 and 1.2 m (Figure 4c) and the optimal size was 0.8 m for this type of landscape (Table 1). By decreasing the value of spikes, small non-ground objects were removed from the final point cloud.

Although the TE obtained at both sites was high because the referenced data included the ground elevation data collected in May, the MF and PTD filters performed satisfactorily compared to TIN (

Table 2. Comparison of the ratio of incorrectly classified points to the total number of points tested (TE) of the filtering algorithms using their optimal set of parameters for the two subscenes, Plot 1 and 2, monitored in June, July, and August.

Subscene/Date	PTD	TIN	MF
Plot 1/26 June	15.57	26.38	9.27
Plot 1/14 July	28.75	46.10	18.28
Plot 1/12 August	23.84	40.56	16.36
Plot 2/26 June	14.75	34.73	19.62
Plot 2/14 July	18.32	37.46	25.65
Plot 2/12 August	16.84	32.39	20.52
Mean Error	19.67	36.27	18.28

). The MF achieved the minimum TE in Plot 1, which was a homogeneous area, whereas PTD performed better in Plot 2, which consisted of low and high vegetation. The nature of the dense vegetation area resulted in a small number of ground points and, consequently, in fewer ground seeds for the interpolation-based filters (TIN and PTD), affecting the densifying process where the morphological filter could provide more ground seed points in general, enabling better coverage. The PTD probably performed better than TIN because TIN considers the lowest point in each grid with a fixed size as ground seed points [48].

3.2. Validation of Canopy Height Models

To validate the resulting CHMs, the field-based measurements of tree height, a building, and plants within randomly selected and geolocated experimental plots (Figure 5a) were compared to the height of the respective objects within each CHM. The heights of these roughness objects were calculated after: (i) segmenting each point cloud dataset into terrain and non-terrain points using the optimized parameters of the MF and PTD segmentation approaches, (ii) normalizing each dataset to the ground (e.g., Figure 5b–d), and (iii) rasterizing the normalized point clouds to generate the CHMs.

Table 3. Comparison between LiDAR-derived height of plants (h) and plants’ h measured manually in geolocated experimental plots. Number of plant samples = 21.

Date	Field h (m)	LiDAR h (m)
26 June	0.52	0.44
14 July	0.71	0.56
12 August	0.78	0.72

outlines the average height of vegetation determined from the field-based measurements in each month and from GIS analysis of the respective LiDAR-derived height values. Regression analysis between the LiDAR canopy height (h_LiDAR) and measured canopy height for 21 plants (h_field) exhibited a linear relationship (Equation (7)) with a coefficient of determination (R²) of 0.89 and root mean square error (RMSE) of 0.028 m.

$${\mathrm{h}}_{\mathrm{LiDAR}}=0.98 {\mathrm{h}}_{\mathrm{Field}}-0.07$$

(7)

The vegetation dynamics in terms of height metrics and the volume of the plants within a subscene area covering 850 m² around the EC tower were calculated for each UAV survey to quantitatively assess the validity of the CHMs (

Table 4. Vegetation dynamics as expressed by the statistics of height of plants (h) and the net volume of a subscene of the CHMs around the EC tower covering 850 m².

Date	Mean h (m)	Mode h (m)	Standard Deviation h (m)	Increase in Vegetation (%)
26 June	0.61	0.60	0.11
14 July	0.76	0.75	0.16	30.25
12 August	0.91	1.00	0.15	44.36

). The estimated CHMs indicated an increase in vegetation height and volume from June to August, but the standard deviation of vegetation height decreased in August. This pattern could be explained by the mechanical removal of potato vines and the ridging of soil to cover growing tubers that both occurred at the end of July to facilitate the harvest of the potato plants by the end of August.

3.3. Source Turbulent Areas Using Morphometric Models

The footprint model of Kljun et al. [66] is a parameterization of a Lagrangian stochastic particle dispersion model and was applied to calculate the extent of turbulent source areas for each 30-min period of the EC observations (Figure 6). For the comparison analysis, we used the morphometric Z₀ values as calculated for the cross-sections of the CHM that coincided with the wind direction that was considered as the input for each run of the footprint model. The footprint model requires the standard deviation of the lateral velocity component, the measurement height, the Obukhov length, the friction velocity, and wind direction (all derived by the EC system), an estimation of the boundary-layer height, and a minimum fetch around the EC tower (approximately 100 m). The 80% cumulative source area for each 30-min EC measurement was utilized to weight the fractional contribution of each grid square of the CHM. This allowed the calculation of a single value of Z₀_EC (by weighting the values in the source area) and the calculation of the average morphometric-derived Z₀ for each turbulent source area (Figure 7). Unstable atmospheric conditions were defined as those corresponding to the ratio z/L < −0.032, while near-neutral atmospheric conditions reflected the relation −0.032 ≤ z/L ≤ 0.032 [67]. The source area climatology was biased toward the dominant west-southerly wind direction, as in Plot 1. During the experimental campaign, only in June did wind originate from the east-southerly direction, as in Plot 2.

3.4. Comparison of Methods to Derive Roughness Length

The mean values of Z₀ as calculated by the anemometric and all morphometric methods gradually increased from June to August, following the progression of vegetation growth resulting from the increases in fai and Z_d due to the higher vegetation density and height. On average, the morphometric-based Z₀ presented strong linear correlations with Z₀_EC (Figure 8), and a standard deviation of less than 4.2 cm with averages ranging from an underestimation of 1.3 cm (Z₀_RT) to an overestimation of 1.9 cm (Z₀_MR) (

Table 5. Comparison of Z₀ (m) derived by the morphometric methods with the anemometric Z₀ averaged over the source areas for daytime hours (from 9:00 to 19:00). The data were screened for wind speed higher than 2m/s and friction velocity higher than 0.2 m/s.

	Z0_RAP	Z0_RT	Z0_MR	Z0_EC	fai
Differences to Z0_EC
June (n = 63)	0.009	0.037	−0.025	0.148	0.048
July (n = 62)	0.018	0.008	−0.023	0.171	0.048
August (n = 36)	0.015	0.013	−0.006	0.200	0.058
Average	0.014	0.013	−0.019
Standard deviation	0.031	0.042	0.022
Unstable conditions Plot 1
June (n = 8)	0.028	−0.003	−0.043	0.117	0.039
July (n = 46)	−0.016	−0.026	−0.042	0.137	0.045
August (n = 3)	−0.022	−0.065	0.022	0.123	0.042
Neutral conditions Plot 1
June (n = 13)	0.078	0.05	0.006	0.171	0.041
July (n = 16)	0.027	0.018	−0.053	0.182	0.042
August (n = 33)	0.018	0.02	−0.009	0.207	0.060
Unstable conditions Plot 2
June (n = 36)	−0.020	0.036	−0.035	0.141	0.054
Neutral conditions Plot 2
June (n = 6)	0.017	0.077	−0.002	0.188	0.046

).

The observed positive correlation between the mean Z₀ obtained by the EC method and the mean Z₀_RT (i.e., the height of the plants) during the vegetation-growing period indicated that an accurate representation of vegetation height derived by a LiDAR system could be effective for estimating Z₀ using the simple rule of thumb method (Figure 8). However, the correlations between Z₀_RAP and Z₀_MR with Z₀_EC exhibited higher coefficients of determination (R² = 0.96 and 0.93, respectively) and smaller RMSEs compared to Z₀_RT. Thus, the investigation of a suitable morphometric method may be crucial to improving the accuracy of canopy aerodynamic characteristics estimations.

The overall difference in Z₀ derived by RAP and the anemometric method was less than 10% for June, July, and August. The RT method had a similar performance to RAP for July and August, while the estimated Z₀_MR was 4% to 19% greater than the mean Z₀_EC. The mean roughness length values under near-neutral conditions were higher than the Z₀ calculated for unstable conditions (Table 5), since the extent of the turbulent source areas was typically larger in the former case with smaller values of friction velocity or wind speed. Perhaps the inclusion of the effect of frontal surface U and u_* in the RAP method as well as the inclusion of vegetation height variability in the MR method enabled the capture of Z₀ amplification under near-neutral conditions, whereas the dependency of Z₀_RT to the averaged h per grid cell produced similar roughness lengths for unstable and near-neutral conditions. The mean Z₀_EC, Z₀_RAP, and Z₀_MR in August and under unstable atmospheric conditions were smaller than the respective Z₀ in July. This could be attributed to the decreased standard deviation of vegetation height within the cumulative source areas observed in August ($\sigma h$ = 0.1 m), which may have translated to a higher density in foliage compared with the respective one for July ($\sigma h$ = 0.15 m).

The role of vegetation structure in the roughness length was identified in the Z₀ variations retrieved in June. The average Z₀ in June calculated for the source areas corresponding to Plot 1 (mean h = 0.61 m) was smaller than that obtained in Plot 2 (mean h = 0.55 m) for both unstable and near-neutral conditions. The mean friction velocity in the direction of Plot 1 was also considerably higher, indicating turbulent disturbance (u_* = 0.53–0.54 m/s in Plot 1 vs. u_* = 0.38 m/s in Plot 2). The decreased Z₀ in Plot 1, even if the average plants’ h was shorter than the respective one in Plot 2, could be attributed to the concurrent lower roughness vegetation density (fai) of Plot 1 and the higher planar vegetation density (pai) compared to Plot 2 (Figure 9), which was the experimental field covered by short grass. The roughness density (fai) is related to the shapes of plant crowns and to the average density of the canopy elements (pai). The pai is related to the effect of intervening spaces between roughness elements in the overall drag efficiency of a canopy, where a higher pai has a smothering effect on the canopy that increases the Z_d; therefore, Z₀ would decrease for a given value of fai, as in Plot 1.

Traditionally, it has been observed that as fai increases, the magnitude of Z₀/$\overline{h}$ also increases at some intermediate level of vegetation densities until it reaches a maximum value for a critical value of fai. The critical fai depends on the method used to determine Z₀ and can be interpreted as the level of homogeneity of the canopy at which adding further roughness elements to the surface does not affect the bulk drag because additional elements merely shelter one another [68]. Similarly, as the pai approximates unity, the surface elements are so densely packed that they merge to form a new surface with limited resistance to airflow. In this study, the spatial distribution of Z₀ exhibited maximum values for pixels with an fai of 0.08 (Figure 10), where the height of the momentum sink starts to move upward since a large fraction of the total drag is exerted by the outermost leaves and branches rather than the background. The drag Z₀ values for pixels with an fai higher than 0.08 are expected to be smaller than the maximum Z₀. This pattern, however, cannot be described with the application of the RT method.

3.5. Influence of Wind Orientation for Deriving Roughness Length

Based on morphometric-based methods, Z₀ can be mapped for a larger area beyond the turbulent source areas of the EC tower by calculating the CHM and the wind directions for which the morphometric parameters are modeled. By considering an isotropic surface, the LiDAR-derived height metrics for each wind direction can be integrated into one value in each grid of the map [47]. Figure 11 illustrates the maps of Z₀_RT, Z₀_RAP, and Z₀_MR for a subscene corresponding to the wind regime of 190° to 347° (130 × 250 m), which was surveyed in June. The RAP and MR tended to exhibit higher roughness length values than the RT. The averaged Z₀_RAP over the homogeneous part of the crop field (e.g., Figure 11) resulted in similar values when different win directions were considered. However, when the area was heterogeneous, the spatial distribution of Z₀ varied significantly.

To assess the influence of wind direction on the spatial distribution of aerodynamic roughness length over a heterogeneous area with low and high vegetation, Z₀_RAP was mapped for a subscene (200 × 300 m) corresponding to the wind regime spanning from 90° to 190° (Figure 12). The orientation of the trees and lower vegetation compared to the wind direction altered the spatial distribution of Z₀, with higher values of Z₀ occurring when the tree arrays and the longitudinal dimension of tramlines were perpendicular to the wind flow (105°). Lower Z₀ values were observed when crops were perpendicular to the wind flow (205°), reflecting the smaller frontal surface of the roughness objects opposed to the wind orientation compared to the frontal surfaces opposed to the wind direction of 105°.

4. Discussion

Morphometric approaches may facilitate the estimation of the aerodynamic surface roughness of an agricultural field with trees at high spatial resolution and over a larger area than the EC turbulent source area but require precise measurements of the height of roughness objects, such those obtained by a UAV-LiDAR system.

The comparison between the height of potato plants and trees as derived by UAV-LiDAR with the respective field-based measurements indicated that UAV-LiDAR may yield strong predictions of the actual height of vegetation canopy. The observed low RMSE is consistent with other reported accuracies for agricultural fields populated by wheat, triticale, and rice [69,70,71] using a terrestrial LiDAR configuration. Within low vegetation ecosystems, airborne LiDARs with a small-diameter footprint exhibited good performance in deriving correlations between the LiDAR-determined vegetation heights and those measured in the field moderately well [72]; however, it was also well-documented that ALS tends to underestimate the vegetation height by at least 30%, depending on the point density [73]. UAV-LiDAR systems may decrease this uncertainty because they generate point clouds with much higher point density, increasing the likelihood for the laser pulses to reach the underlying terrain and reflect from the maximum top of every stem they contact. However, the proper classification of ground and vegetation point cloud data is a crucial step for further processing the geometric characteristics of a canopy in order to calculate the aerodynamic and biophysical properties of low-vegetation canopies. For example, Wang et al. [74], using a slope- and angle-based filtering method [75] to classify UAV-LiDAR point clouds into ground and grassland, found that LiDAR underestimated the canopy height, whereas at the locations where the grassland was less than 5 cm, LiDAR-derived heights were overestimated. To separate ground points from points representing winter wheat as acquired by UAV-LiDAR [33], the authors argued that the cloth simulation filtering algorithm [76] had to be parameterized according to the temporal variations in the vegetation density. Since criteria for choosing the most appropriate filtering method and the optimized parameters of every filtering algorithm are lacking, we used the calculation of the total errors of morphological and interpolated-based filtering algorithms by setting different values of their parameters. This approach eliminated misclassification and the resulting biases in the CHMs without requiring the visual interpretation of the results, which could consider the effect of heterogeneous vegetation on the segmentation process for this specific land cover and topography. The morphological filter could detect more ground points in the dense low vegetation site, whereas the PTD performed better when trees and low vegetation had to be segmented. In both cases, the optimal window size needed to be considered, which depends on the sizes of the contained objects [77] and is subjected to variations in the acquired point density, which can be altered by changing the flight line distance or the flight speed.

Another critical issue for defining Z₀ at both high temporal and spatial resolution using UAV-based observations is the appropriateness of the applied morphometric model. In agricultural studies that are based on remote sensing data to estimate turbulent heat and gas fluxes, Z₀ is usually calculated as a fraction of the mean height of roughness objects (e.g., [78]) or as a dependent variable of a VI assuming a homogeneous area to derive the resistance to heat transfer [4,5,6,7]. However, the temporal variations in Z₀ even in agricultural areas do not always simply correlate with the vegetation height, or the area is not always homogeneous; therefore, these approaches may introduce uncertainty to the estimated energy or gas fluxes [79]. In this study, Z₀ depended on the height-based geometric parameters and on the vegetation structure, expressed by the frontal and planar area indices (e.g., Figure 9 and Figure 10), resulting in better correlations between Z₀_RAP and the Z₀_EC averaged within each source area compared to Z₀_RT (Table 5).

The mean Z₀_MR could be considered as the upper limit of the mean Z₀_EC using a 0.25 m subcell scale that can describe the height variation of a plant. This method can accurately capture the vegetation variability using high-resolution CHMs (e.g., 0.10 m) from UAV-LiDAR measurements, but it may be sensitive to the choice of the grid cell size and its subcell sizes. For example, it was observed that coarser subcell scales resulting from coarser ALS-derived CHMs can lead to less convergence between Z₀_EC and Z₀_MR [25]. The RAP method may be more appropriate to simulate the temporal variation in Z₀ in this type of landscape, since the formula describes an interplay between roughness length and alterations in airflow orientation (through fai) and wind speed for regular arrays of roughness elements that characterize our study site.

The assessment of the accuracy of the Z₀_RAP and Z₀_MR values across different wind fields outside the turbulent source areas would necessitate the acquisition of anemometric-derived Z₀ across the whole field. Therefore, a precise statement about how these morphometric-based Z₀ respond to vegetation variability is difficult. However, differences in the spatial distribution of Z₀ between crop fields, bare soil, and trees were generated (e.g., Figure 11), and the effect of wind direction on the fai and Z₀ for the heterogeneous subscene of the agricultural site was evident using the RAP method (Figure 12). These observations are aligned with the findings of Colin and Faivere [23], who documented that the RAP approach could account for the heterogeneity of an area covered by grassland with staggered arrays of trees.

The morphometric models, though, do not account directly for the potential effect of vegetation’s porosity on the amount of drag exerted on the flow [80]. Thus, it could be claimed that the calculated aerodynamic resistance would be overestimated by considering the porosity of the foliage structures equal to one. That would probably be more effective for areas where the dominant roughness element is trees, which, in contrast to compact obstacles, mainly oppose a resistance to the airflow for tree heights with the highest foliage density. Kent et al. [81] found that the effect of the porosity of low vegetation on the roughness length compared to higher obstacles was negligible. For more complex croplands, a sensitivity analysis of the drag coefficient used in various morphometric models may provide further information regarding the alterations in shear stress sourced from the internal structure of high vegetation. Although we tested the applicability of the Macdonald [12] and Millward-Hopkins [82] morphometric models, which can account for the differential drag imposed by different types of obstacles, the results were not exploited in this study because they did not capture Z₀ variations in low vegetation, probably because these models were designed to determine Z₀ in urban areas.

The methods assessed here could generate the spatial distribution of Z₀ at centimeter-level resolution for an agricultural site and for selected prevailing wind directions, but the contribution of the upstream roughness elements cannot be quantified. The choice of an appropriate spatial scale of analysis for computing Z₀ could be derived from a general analysis of the landscape characteristics along the prevailing airflows.

5. Conclusions

In this study, we explored a method of generating maps of roughness length at ultra-high spatial resolution (0.10 m) of a typical dense temperate agricultural field with sparse trees using the newly developed UAV-LiDAR scanners and morphometric roughness models without requiring any other data sources, such as optical photogrammetry or eddy covariance observations. This method is particularly useful for enhancing the sustainability of farming practices and recreation activities in agroforestry and agroecology applications, since the spatial distribution of Z₀ can delineate how land cover alterations affect shear stress and turbulence, which, in turn, regulate the air–surface exchange of energy, water, and greenhouse gases.

For the determination of Z₀, the selection of the appropriate morphometric method and the effective classification of UAV-LiDAR-derived point clouds into vegetation and terrain that generates precise CHMs are both critical.

Overall, convergence was observed between the EC-derived Z₀ values and those found through the UAV-LiDAR-driven models at the turbulent source area scale. All morphometric models showed a standard deviation of less than 4.2 cm with averages ranging from an underestimation of 1.3 cm (Z₀_RT) to an overestimation of 1.9 cm (Z₀_MR). The detailed comparison indicated that the Raupach roughness model is more suitable for simulating the temporal variations in Z₀. The spatial distribution of zo_RAP for a heterogeneous subscene beyond the turbulent source areas was conditioned by the shape of the frontal surface opposed to wind direction, with a higher Z₀ occurring when the tree arrays were perpendicular to the wind flow.

A sensitivity analysis of three filtering approaches to segment the point cloud data to low vegetation and ground highlighted the errors associated with CHM preparation. The morphological filter performed satisfactorily over the more homogeneous area covered by plants with a 1 m window size of the opening operation while the filter was not so sensitive to elevation threshold due to the relative flat landscape. The PTD filter produced fewer errors in a subscene consisting of low and high vegetation for an iterative distance close to 0.4 m and an iterative angle of 4°. The TIN interpolation-based filter generated more errors in detecting ground and non-ground points for this type of vegetation and landscape.

This technique may advance the establishment of more precise spatial representation of potentially nonlinear relationships between canopy structural characteristics and surface water and gas dynamics across heterogeneous land cover types. The application of all the elaborated approaches to derive precise CHMs and Z₀ values in agricultural fields with other crop species and different climatic conditions could be used to assess their adequacy in other contexts. Further research is needed to improve the morphometric models for Z₀ in vegetated landscapes that can benefit from canopy height models of ultra-high spatial resolution to account for surface drag effects of upstream roughness elements.