Quality of digital elevation models obtained from unmanned aerial vehicles for precision viticulture

Methods and results: Experiments were carried out on a 4-ha vineyard located 10 km north of Beziers (France). The experimental site presents slope and aspect variations representative of mechanised commercial vineyards in Languedoc Roussillon. DSMs were provided by three UAV companies selected for the diversity of their solutions in terms of image capture altitude, type of UAV and image processing software. DSMs were obtained by photogrammetry and correspond to commercial products usually delivered by UAV companies. DSMs from UAV were compared to a reference Digital Elevation Model (DEM) acquired by a laser tachymeter. Four indicators were used to test the quality of DSMs: the mean error and its dispersion in the XY plane and in elevation Z. Results show a good georeferencing of the DSMs (MeanErrorXY<10 cm) and a similar quality in elevation (MeanErrorZ<10 cm) estimation. Results also show that the error in elevation is highly spatially structured. The spatial patterns observed did not depend on the elevation and could be related to algorithms used to compute the DSMs.


INTRODUCTION
The detailed knowledge of elevation and its variations is of importance in viticulture.These data are typically contained in information layers called Digital Elevation Model (DEM).DEMs with high spatial resolution (<1 m) and high accuracy (<0.05 m) in elevation measurements are interesting sources in Precision Viticulture (PV).Indeed, several studies have shown the relevance of this information in understanding and explaining the spatial variability of the production parameters like yield, vigour and vine water status, either at the whole vineyard scale (Santesteban et al., 2013) or at the within field scale (Bramley et al., 2011a).Beyond elevation, DEMs also make it possible to estimate the slope and the aspect of the fields.Therefore, DEMs provide basic information for soil erosion investigations (Martinez-Casasnovas and Sánchez-Bosch, 2000), delineation of optimal maturation zones (Olsen et al., 2011) and terroir (Carey et al., 2008), as well as identification of frost hazard risk areas (Madelin and Beltrando, 2005).Until now, the implementation of an accurate DEM (precision<0.05m) with a high spatial resolution required conducting a field survey with specific sensing systems either manually (Santesteban et al., 2013) or embedded on mobile machines (Bramley et al., 2011b).Photogrammetry by Unmanned Aerial Vehicle (UAV) (Stefanik et al., 2011) is a new source of information that allows the production of Digital Surface Models (DSMs) with a high spatial resolution.This technology may present an alternative to conventional methods.DSMs provided by commercial services based on UAV are already available for the wine industry.However, the use of this source of information still raises some questions.Indeed, UAVs do not necessarily measure the elevation of the ground, but the elevation of the visible objects from above (i.e.vegetation over the soil).This justifies the distinction between DSM and DEM.
The interest of UAV-based DSM has been shown for applications in different areas like archeology (Sauerbier and Eisenbeiss, 2010), hydrology (Leitão et al., 2016), engineering (Uysal et al., 2015) or agriculture to monitor crops (Bendig et al., 2013).These studies showed it was possible to estimate the elevation with errors ranging from 0.5 m to 0.05 m.However, they were all performed with specific equipment and acquisition conditions as part of a scientific work.The assessment of the quality of DSMs is made point-to-point on a small amount of ground truth sites (Uysal et al., 2015) or by type of observed object (Leitão et al., 2016) on areas where characteristics (slope, orientation) and objects may differ drastically from that of a vineyard context.
Considering the requirements of the wine industry, the objective of this work is to assess the quality of elevation measurements derived from DSMs obtained by UAVs in the specific context of a Mediterranean vineyard.Given the wide variety of parameters (type of UAV, acquisition conditions, chip size of the sensor, flight elevation and speed of the flight, image overlap and software and algorithms used) that may affect the quality of the results, this work does not aim to provide references on the best acquisition conditions.It focuses on comparing several current commercial services making the assumption that the companies in charge of these services deliver the best possible information under the conditions of the study.
The originality of this work is: i) to consider several acquisition and processing chains which are currently available for DSM by a UAV: type of UAV, elevation of image acquisition, software used to compute the elevation from images, ii) to consider a study site as representative as possible of vineyard landscapes; this site takes into account specific objects which need to be described by the DEM, it also shows the specific constraints in obtaining a DEM from a DSM in these specific conditions, and iii) to propose a detailed study of the spatial distribution of the error over the whole study area in order to verify whether the error is constant over the entire study site or whether specific objects of viticulture landscape affect the quality of the DEM derived from UAV.To our knowledge, such a study has never been carried out in a commercial context.

Study zone
The study zone was a 4-ha vineyard located 10 km north-west of Beziers, in Languedoc in the south of France.It was chosen for its representativeness of the mechanised commercial Languedoc vineyards.The four fields were around 1 ha each.The elevation ranged from 150 m to 175 m.As shown in Figure 1, the study zone was organised as a basin with an outlet on the south-east corner.It presented different slope orientations.The main slope ranged between 0 and 10% with a north/south orientation.The characteristic features of the landscape were a zone of break in slope with an embankment running from east to west, and two paths above and below the bank, corresponding to flat zones.The four fields were non-irrigated, established with 1-m spacing between vines and 2.5-m inter-rows and trained in a vertical shoot positioning system with no cover crop (mechanical weeding).

Reference measurements
To enable georeferencing of all the information measured on the study zone, seven reference points were defined on the study zone (Figure 1).These reference points were positioned on visible targets with position and location measured using a RTK GPS receiver (error<0.02m) (Topcon Hiper SR, Livermore, United States).
Elevation was measured accurately on 390 sampling points over the study zone (Figure 1).This sampling scheme was carried out following a grid of approximately 20*20 m over the study zone.This grid was supplemented with measurements at specific locations corresponding to particular changes in topometry like hydrographic network or embankments.At each sampling point, elevation was measured using a laser tachymeter (TPS 1001, Leica, Heerbrugg, Switzerland) with a measurement error less than 0.01 m.A DEM (Figure 2b) was then derived by Kriging elevation values on a 1*1 m grid.The resulting resolution of the DEM derived from the tachymeter (considered as the DEM of reference) was therefore of 1 m.
Among the 390 sampling points, 48 were positioned on bare soil next to objects like trellising posts or irrigation boxes, easily identified on UAV images.These points were used to match the reference coordinates and the coordinates of images from UAV, allowing estimating the quality of georeferencing.In order to be able to compare the time required for obtaining the DEM of reference with those required by UAVs, the acquisition points throughout the study area needed 2 people for 7 hours (14 hours) and 2 hours for data processing (kriging and mapping).

Data acquisition
DSMs were provided by three commercial companies: A, B and C.They were obtained by photogrammetry on images acquired by UAV with fully commercial data acquisition chains.The corresponding orthophotos in the visible range (RGB) at the same resolution were also supplied.
As shown in Table 1, the three acquisition chains were based on the two types of UAV used in commercial applications: multirotor (company A) and flying wing (companies B and C).
Flight altitude and distance to pilot were determined according to the scenario of the French legislation, which defines operational scenario for commercial use of UAV according to the maximum flight elevation, the maximum distance to the pilot and the maximum embedded weight (Journal Officiel de la République Française, 2015).Two scenarios were used: S1 scenario for companies B and C and S2 scenario for company A. Images were georeferenced from either the reference points presented above (companies A and C) or from points measured by the UAV company itself (company B).Post-treatment was made with softwares supplied as standard by UAV constructors (Table 1).Spatial resolution of images corresponded to the commercial products usually provided to customers by each company.

Comparative method
To characterise the accuracy of the three DSMs, they were compared to the reference DEM on 48 comparison points.The comparison points were identified on each of the three corresponding orthophotos by photo interpretation using the GIS software QGIS (Quantum GIS Development, V2.8, Open Source Geospatial Foundation).Note that this comparison approach aimed at characterising absolute error of DSM.It assumed that if errors exceed a given threshold specified by the case study under consideration, then elevation data from UAVs could hardly be used to produce other derived information such as slope or orientation.
For each DSM and each comparison point, the X and Y coordinates were extracted from the corresponding orthophoto and the Z coordinates were extracted from the DSM.The X, Y and Z values of the 48 comparison points where then compared for each DSM to the reference DEM.
Four descriptive indexes were used to compare the three DSMs to the reference DEM: MeanErrorZ, MeanErrorXY, RMSE_Z, and RMSE_XY.The MeanErrorZ and the MeanErrorXY estimated, in average, how close to the DEM of reference are the DSM in Z and in the XY plane respectively.RMSE_Z and RMSE_XY estimated the deviations of DSM when compared to DEM of reference in Z and in the XY plane values respectively.Indexes were described in Equations 1, 2 and 3.
Where: j refers to the DSM (with j = 1:3 ) i refers to a given reference point (with i = 1:48) n refers to the total number of reference points (with n = 48) Z DSM i (j) corresponds to the estimated elevation at a given point from a given DSM Z ref i corresponds to the reference elevation at a given point [Eq. 2] MeanErrorXY and RMSE_XY are computed by replacing the term ErrZ i (i) with ErrXY i (j) in Equations 1 and 2 respectively.ErrXY i (j) is presented in Equation 3.
-104 - OENO One, 2016, 50, 3, 101-111  ©Université de Bordeaux (Bordeaux, France) Léo PICHON et al. Where: X DSM i and Y DSM i refer to the coordinate of a given point from a given DSM X ref i and Y ref i and refer to the coordinate of a given point on the reference DEM

Expert classes considered for the errors
To highlight the potential applications of DSM in viticulture, different classes of errors were considered.These classes were defined from expert specifications developed by a company specialised in agri-environmental studies in agriculture and viticulture (Envilys, Villeneuvelès-Maguelone, France; Delalande, personal communication).In this study, these classes were used to define possible uses for DSM obtained from UAV.These classes were indiscriminately applied to the average error (MeanErrorZ) and punctual errors (ErrorZ); they are detailed (in absolute value) hereafter: >0.2 m: inadequate for applications at the field level, <0.2 m: suitable for slope computation or field aspect estimation, <0.1 m: suitable for embankment or ditch identification and measurement, <0.01 m: suitable for within field soil erosion or gully characterisation.
Considering the accuracy of our reference measurement (0.01 m), this last class will only be considered to highlight the error detection limit of our method.
The requirements on ErrorXY and MeanErrorXY depend on the need to retrieve one exact place on the field or to overlap other information layers.These issues were not considered as specific to this study so no specific classes were used.

Interpolation and mapping
Data mapping and photointerpretation were performed using QGIS (Quantum GIS Development, V1.8, Open Source Geospatial Foundation Project) by importing eastern and northern values for each elevation variable.Data interpolation was performed using GeoFIS (V0.1, INRA/IRSTEA/Montpellier SupAgro, France).The interpolation method used in this study was based on Ordinary Kriging.Semi-variogram analysis and interpolation were conducted with GeoFIS.ErrorZ maps were computed following the same method.

RESULTS
On a first approach, at the study zone scale, the three DSMs are very similar to each other.The main landscape components are clearly visible (Figure 2a): the bank running from east to west, the north southern slope and the basin with an outlet on the south-east corner of the study zone.Those elements are also visible on the reference DEM (Figure 2b).On the edges, some areas are not covered by the reference DEM due to the interpolation method.Repeatability tests (results not shown) result in an error in the reference measurements lower than 0.01 m.This result is in agreement with the accuracy announced by the tachymeter manufacturer.It confirms the limits of our experimental setup to characterise errors less -105 -OENO One, 2016, 50, 3, 101-111 ©Université de Bordeaux (Bordeaux, France) When focusing on a particular area (Figure 3a) with vegetation and a break of slope (Figure 3b), dissimilarities between DSM and DEM appear.
The vegetation is clearly visible on the DSM (Figure 3c) and one can even almost perceive the shape of the canopy.This fact highlights the potential benefits of investigating how information should be extracted from DEMs from UAVs in order to improve their performance in PV.
Figure 3d shows triangle patterns along the bank on the reference DEM.Those patterns are artefacts caused by the locations of the measurement points and the interpolation method.Concerning the elevation, similar results were obtained (Figure 5a and Figure 5b).Indeed, all DSMs provided a really good and identical deviation (<0.1 m), but a significant shift on the  error for company B (>0.7 m) compared to companies A and C (<0.1 m) was observed.The shifts observed either in Z or in XY for company B may be due to the choice of the reference points used to georeference the information.Indeed the company used its own reference points instead of those provided within the framework of study.Note, however, that this choice does not impact the dispersion (RMSE) of the data either in Z or in XY when the shift is removed from the data.
In order to analyse the spatial distribution of errors (ErrorZ, Eq. 1), an analysis of the spatial autocorrelation was performed with a semi-variogram model.Table 2 shows the parameters of the semivariogram models of the errors for the three DSMs.The nugget effect (C0) corresponds to an erratic phenomenon, the Sill corresponds to the total variance of the error and the range indicates the distance over which the error is no more spatially auto-correlated (i.e.beyond this value, error values become independent).The difference between Sill and C0 is the variance corresponding to a spatially organised phenomenon (also called C1).
The three DSMs present the same C0 value (Table 2).It corresponds to a standard deviation of 0.03 m.This result shows that an erratic phenomenon affects the quality of the elevation values with the same magnitude for the three DEMs.This erratic phenomenon may result either from error due to the measurements itself (data acquisition, data processing) or from the quality of the reference measurements used in this experimentation.This point will be discussed later.
For the three DSMs, a significant proportion of the variance corresponds to a spatially organised phenomenon.It accounts for almost 50% of the total variance for DSM B and DSM C and for more than 80% for DSM A. Moreover the range shows that for the three DSMs, this spatial phenomenon occurs over a large distance (>150 m).This shows that a significant part of the error does not correspond to a random phenomenon but to a spatially organised phenomenon.This observation points to the presence of large spatial patterns where errors are higher or lower on the study area.This results in large areas where elevation data are systematically over or underestimated.This phenomenon may lead to the presence of artificial gradients which can be problematic to estimate the slopes over distances greater than 150 m.In order to better analyse this spatial organisation and to verify any relationship with specific characteristics of the study area, error maps were performed (Figure 6).
Figure 6 confirms the results of the previous section.Clear spatial patterns are observed whatever the DSM considered.Note that spatial patterns are different from one DSM to another showing that the origin of this phenomenon is not due to a systematic problem in reference measurement acquisition.Figure 6 shows clearly that error values are not randomly distributed over the study zone.
The lower vertical accuracy previously observed with the DSM provided by company A (Figure 6a) is explained by two zones on the western and eastern parts of the study area where error values are higher than 0.2 m (in absolute value).For the other DSMs (Figures 6b and c), no value higher than 0.2 m is observed.This result shows that in our experimental conditions, DSM B and C may well be used for slope computation and field aspect determination.Both these DSMs may be suitable for embankment identification and ditch measurement but on a smaller area corresponding to error values <0.1 m.For DSM A, only the central part of the study zone presents the required quality for the considered applications.
This experiment does not allow the cause of spatial pattern errors to be identified.Indeed no relationship between errors and elevation values (results not shown) was observed.Therefore higher errors may occur either in the low or high parts of the study zone.The spatial organisation of errors may be related to the 3D construction algorithm and its implementation.This known phenomenon (Schleicher et al., 2007) corresponds to an accumulation of errors on image sequences during image processing.It involves a complex phenomenon which depends on the flight route, image overlap and the selection of reference points.This could for example explain the higher errors observed over the entire eastern part of DSM A.
From a practical standpoint, this result shows that assessing the quality of DSM from UAVs must always be based on reference points distributed over the whole study zone under consideration.It also shows that errors are independent from the shape or elevation variations of the study zone, but are spatially auto-correlated

DISCUSSION
With current commercial solutions, photogrammetry by UAV may provide accurate elevation measurements that can answer the demand for relevant applications for the wine industry.Indeed, under the conditions of the experiment, DSM from UAV complies with significant applications for terroir characterisation and vine field classification.In particular, it is possible to estimate the elevation of the fields precisely with a resolution relevant enough for slope and aspect estimation.Therefore, UAV may constitute a relevant sensing platform for these applications providing a more affordable solution than conventional methods which require a specific ground survey (RTK GPS, laser tachymeter) performed by a specialist.It also constitutes a more reliable source of information than most available DEMs which can be free such as Aster MNT (resolution 30 m, ErrorZ ~ 2 m.) or not free such as IGN-RGE ALTI (ErrorZ ~ 0.5 m) (Devaux, 2015).
In our conditions, the accuracy achieved is however lower than the accuracy of the reference measurement (0.01 m) on a significant part of the vineyard.As already mentioned, problems with 3D construction algorithm are certainly the main cause of errors.However, three other sources of inaccuracy have to be considered: -The conditions of image acquisition are not optimal enough.Many factors may affect the quality of the DSM provided by UAV (flight altitude, chip size of the sensor, image overlap during acquisition, software and algorithm used to find the corresponding points between images).
Investigating the incidence of all these parameters together was not the purpose of this study, which assumed that commercial companies mastered their acquisition chain to provide the best possible information.It is possible that a flight at lower altitudes (and thus with a higher spatial resolution) would result in a DSM with better accuracy.Note however that these flight characteristics were carried out with company A (Table 1) and did not achieve better accuracy in our conditions.Moreover, this high resolution acquisition necessarily leads to a significant volume of data.
Assuming a lower flight altitude may improve the quality of elevation measurements, these conditions entail greater constraints that seem currently difficult to manage for a commercial vineyard.These conditions could however be considered in the context of research projects with a partner specialised in the processing and storage of large volumes of data.
-The error in XY which may affect the spatial correspondence between the reference measurements and the Z values of the DSM from UAV.These differences in XY may have led to compare Z values that are not co-located, which may increase the observed error in Z.
-Finally, the quality of the reference measurements used is also an important limitation in our study.Repeatability tests resulted in an error in the reference measurements lower than 0.01 m.However, semi-variogram analysis (Table 2) highlighted an erratic error of 0.03 m.Assuming this erratic error is only due to the reference measurements, two additional sources of inaccuracy may be considered.First, the soil roughness due to mechanical tillage of the interrow may have led to clod formation and therefore might have introduced short range variability in the reference measurements mainly due to the manual positioning of the target of the laser tachymeter.Second, the location of the reference points which was performed manually on the images.This operation is likely to induce an inaccuracy in the positioning of the reference point and a resulting error.Note, however, that if effective, these sources of error result in a random phenomenon.They do not explain the spatial organisation of the observed errors which account at least for 50% of the error distribution.
This study focused on the point-to-point comparison of reference elevation data and data measured by UAVs over a representative study area made of several vine fields.The assumption of this work is that all information (slope, orientation, etc.) calculated at this scale are relevant when precision is met on the whole study area.However, locally and relatively, high resolution elevation information provided by the DSM from UAV seem relevant.This is shown in Figure 4, where breaks in slope areas were better highlighted by the DSM from UAVs than the reference DEM.This observation puts in evidence the potential value of DSM from UAV to identify variations of brutal change in elevation over short distances.Even if the absolute altitude of this phenomenon is imprecise, abrupt changes could be identified and characterised easily.This aspect would require a specific experiment in order to check the consistency of DSM data over short distances.
Regarding PV, DEM derived from UAV image acquisition offers interesting opportunities.Indeed, the spatial resolution of elevation measurements is homogenous with those of other information sources currently used to characterise within field variability.Santesteban et al. (2013) showed the usefulness of relevant elevation data to delineate zones either at the within field level or at the whole vineyard level.Therefore, in association with data layers estimating the vegetative expression either from remotely sensing or proximal sensing techniques or soil characteristic (apparent soil conductivity or resistivity), DEM from UAV may be a convenient additional information source in PV.However, from a practical point of view, the electrical conductivity of soil surveys usually lead to DEM with a similar spatial resolution thanks to an ATV (All-Terrain Vehicle) embedded differential GPS.In this particular case, the realisation of a DSM by UAV does not seem necessary.
The quality of the results obtained allows considering the measurement of the canopy height.This aspect is already being investigated for other crop productions, like olive trees (Torres-Sanchez et al., 2015b) and, more recently, vines (Burgos et al., 2015), reporting quite interesting results when compared to ground reference data.Canopy height from UAVs could provide complementary information to measurements provided by ground sensors like the GreenSeeker® (Drissi et al., 2009).Both measurements could provide estimates of the thickness (with porosity) and the height of the canopy offering the possibility to estimate the exposed leaf area.The information on the canopy height could also be complementary to the index of vegetation whether measured by airborne devices or satellite.Note that the vegetation index can also be measured by UAV, provided it is equipped with a suitable sensor.This feature highlights the value of the UAV platform which, in a single pass or in two successive passes, would simultaneously provide information on the vegetative expression and the canopy height.
Finally, in this work, the DEM was derived from the DSM by a manual identification of sites corresponding to bare soil.Therefore it is assumed that image segmentation of the bare soil and canopy is possible.Such image segmentation does not present major issues.Torres-Sanchez et al.
(2015a) have already shown the possibility to segment bare soil from visible images, provided the resolution of the images is adapted to the image processing.However, in viticulture, applying such image processing techniques assumes that a large proportion of bare soil is visible on the images.This constraint determines the optimal period of image acquisition.Winter or early spring could be the best time of acquisition to obtain a DEM in viticulture.Note, however, that this constraint may be very limiting in the case of winter cover crops.Note also that this characteristic limits drastically the possibility to derive a DEM from DSM in the case of vegetated areas such as embankments or ditches.

CONCLUSION
This work shows that DEMs derived from UAV images are a reliable source of information in viticulture.The DEMs currently provided by UAV companies have an accuracy that can address important issues for the wine industry: for specific problems of precision viticulture to support the delineation of within field zones, for field selection or at a larger scale for terroir characterisation purposes.
However, this work highlights some limitations of this source of information.Indeed, according to the conditions of acquisition, elevation errors are not uniform over the entire study area.Locally, the importance of these errors may limit the use of this information.In particular, for runoff and erosion studies, which require very high accuracy on elevation data, DEM from UAV images may limit, at least locally, the relevance of the results derived from this source of information.
Finally, it is necessary to consider that photogrammetry by UAV images provides digital surface models and not directly a digital elevation model.The estimation of a DEM from a DSM can be difficult in certain conditions when bare soil is hardly visible, especially in the presence of a permanent cover grass or when images are acquired in full vegetation.

Figure 1 -
Figure 1 -Presentation of the studied zone and location of the reference points.

Figure 2 -
Figure 2 -Elevation maps of a) Digital Surface Model (DSM) obtained by UAV (company C) and b) reference Digital Elevation Model (DEM) obtained by the laser tachymeter.

Figure
Figure4aprovides an overview of all results concerning the characterisation of XY position with MeanErrorXY and RMSE_XY of each DSM when compared to the reference DEM.Figure4bsummarises index values.While for the deviation very similar and satisfactory results were obtained for the three companies (<0.3 m), a slight shift was observed on the error for company B (>0.3 m) compared to the other two (<0.10 m), which could be explained by the choice of their own reference points.

Figure 3 -
Figure 3 -Comparison between DSM and DEM on a slope break zone: a) location of the zone, b) orthophoto map of the zone, c) Digital Surface Model (company C), d) reference Digital Elevation Model obtained with the laser tachymeter.

Figure 4
Figure 4 -a) Graphic representation of MeanErrorXY (points) and RMSE_XY (circles) of DSM geopositioning and b) summary table of the corresponding values.

Figure 5
Figure 5 -a) RMSE_Z for the three different DSMs and b) summary table of the indexes used to characterise the elevation error.

Figure 6 -
Figure 6 -Observed error (ErrorZ) spatial distribution for the three different DSMs (Digital Elevation Models): a), b), c) respectively for companies A, B, C.