**Geometric and Radiometric Consistency of Parrot Sequoia Multispectral Imagery for Precision Agriculture Applications**

#### **Marica Franzini 1,\*, Giulia Ronchetti 2, Giovanna Sona <sup>2</sup> and Vittorio Casella <sup>1</sup>**


Received: 29 October 2019; Accepted: 2 December 2019; Published: 5 December 2019

**Abstract:** This paper is about the geometric and radiometric consistency of diverse and overlapping datasets acquired with the Parrot Sequoia camera. The multispectral imagery datasets were acquired above agricultural fields in Northern Italy and radiometric calibration images were taken before each flight. Processing was performed with the Pix4Dmapper suite following a single-block approach: images acquired in different flight missions were processed in as many projects, where different block orientation strategies were adopted and compared. Results were assessed in terms of geometric and radiometric consistency in the overlapping areas. The geometric consistency was evaluated in terms of point cloud distance using iterative closest point (ICP), while the radiometric consistency was analyzed by computing the differences between the reflectance maps and vegetation indices produced according to adopted processing strategies. For normalized difference vegetation index (NDVI), a comparison with Sentinel-2 was also made. This paper will present results obtained for two (out of several) overlapped blocks. The geometric consistency is good (root mean square error (RMSE) in the order of 0.1 m), except for when direct georeferencing is considered. Radiometric consistency instead presents larger problems, especially in some bands and in vegetation indices that have differences above 20%. The comparison with Sentinel-2 products shows a general overestimation of Sequoia data but with similar spatial variations (Pearson's correlation coefficient of about 0.7, *<sup>p</sup>*-value <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16).

**Keywords:** geometric consistency; radiometric consistency; point clouds; ICP; reflectance maps; vegetation indices; Parrot Sequoia

#### **1. Introduction**

#### *1.1. Key Topics*

Precision agriculture (PA) [1] is a very significant societal challenge and promises to enable several significant improvements: increase of productivity; optimal, and thus reduced, use of pesticides and fertilizers; and decreased use of water. These will translate into substantial benefits, including making more food available for mankind, increasing environmental sustainability, and contributing to the mitigation of climate change effects [2]. One key component of precision agriculture is crop health diagnostic capability. Within this context, in the last 5 years the use of lightweight unmanned aerial vehicles (UAVs) equipped with multispectral sensors has become quite popular. UAV-based surveys offer unprecedented ground resolution and operational capability. The second feature is particularly significant when periodic monitoring has to be performed, as the operator is free to choose the optimal time to fly.

#### *1.2. Background*

The processing of large datasets that cover wide areas and which need to be acquired by several UAV missions is still a challenging task. As for photogrammetric projects, these types of datasets, composed of various sub-blocks, require a careful assessment of the accuracy of the final products, from both geometric and radiometric points of view. This is even more true when time-series are analyzed; the consistency between data is mandatory in these cases.

Regarding geometric issues, the recent evolution of UAVs has provided low-cost systems with direct georeferencing (DG) capability. DG has several advantages: it allows flights in remote areas where access could be difficult or impossible [3], and reduces mission time and costs, since no ground control points (GCPs) need to be installed and measured. Unfortunately, navigation-grade GPS/GNSS receivers (Global Positioning System/Global Navigation Satellite System), such as the one integrated with the Parrot Sequoia sensor [4], are not of sufficient quality in the solution position for georeferencing of images. Some authors investigated this topic [3,5,6], obtaining metric errors. GCPs are traditionally suggested for georeferencing purposes. The number and distribution of GCPs have been explored by several authors [7–10]. Independently from the extent of the surveyed area, they state that a small number of GCPs is useful when they are only needed to perform datum transformation, while a larger number is necessary when camera self-calibration must be performed. For this second aim, their spatial distribution is important too, as ground points must cover the whole area of interest. However, GCPs cannot always be guaranteed in some applications, such as in precision agriculture, where inaccessibility is a frequent condition due to crops' stages of growth. In this case, other ground information can be useful, such as pre-existing orthophotos [11].

Independently from the strategy used for images orientation, the geometric quality of the results must be assessed both in terms of accuracy and consistency. Such analysis can be focused on exterior orientation parameters (EOPs) or on the photogrammetric products, such as dense point clouds or orthophotos. EOPs are traditionally evaluated by using a set of check points (CPs), which are considered during the bundle block adjustment as simple tie points; residuals between the photogrammetrically obtained object coordinates of markers and those preliminarily determined by surveying are then evaluated [12–14]. Photogrammetric products are assessed using additional information acquired by alternative systems, such as GNSS receivers and total stations [15,16] or light detection and ranging (LIDAR) [17–19].

Data quality can also be assessed in terms of consistency. This term means the agreement between different (partially overlapping) datasets acquired at the same time or at different times. Consistency can be assessed on various photogrammetric products, such as point clouds or orthophotos. The mentioned criterion is particularly significant when time-series are processed or when different processing strategies are tested, as in our case. Within this framework, [20] evaluated volumetric changes of a landslide areas using point clouds over a time-series, while [21] assessed the consistency of UAV-derived point clouds in relation to the focal length and target set.

While geometry is almost always considered when quality assessment is performed, radiometry is less often investigated but plays a key role in several applications, such as precision agriculture and environmental pollution detection. Regarding radiometry, some critical issues remain unsolved, such as which corrections must be considered and modelled, and several authors recently started to investigate these aspects. Honkavaara et al. [22] studied and assessed a processing methodology for biomass estimation in agriculture with a lightweight UAV spectral camera under varying illumination conditions. In [23], authors captured images from an UAV with Parrot Sequoia and assessed canopy reflectance consistency in avocado and banana orchards in Australia, while in [24] reflectance anisotropy of potato canopies in the Netherlands was mapped with a frame camera mounted on an UAV.

Although in recent years sensor manufacturers have improved in describing sensor performance and providing tools for performing radiometric corrections [25], the radiometric quality of data is still uncertain. The reliability of spectral information acquired by multispectral sensors mounted on UAVs is not completely clear [26]. Absolute accuracy might be insufficient for some applications, so that calibration procedures will be required [27].

The proposed radiometric calibrations are based on the availability of spectral targets, whose reflectance response is measured in situ with a spectrometer [28–33]. After evaluating costs, surveying and processing times, and required instrumentations and expertise, this calibration methodology cannot be adopted as a routine method in precision agriculture. Indeed, farmers must take into consideration data quality as well as economic impact [26]. To ensure the dissemination and the use of calibration procedures in the agriculture sector, it is advisable to optimize the exploitation of consumer-friendly tools, and best practices must be simplified.

#### *1.3. Motivation*

Nowadays equipment vendors are making an effort to supply easy-to-use HW (Hardware) and SW (Software) so that crop monitoring can be performed by individual farmers. The bundle of Parrot Sequoia© (Parrot S.A., Paris, France) and Pix4D© (Pix4D S.A., Prilly, Switzerland) is a clear and popular example of this approach. The present paper arises from one simple yet crucial question: what is the reliability of the radiometric information and of the related vegetation indices acquired by the Sequoia camera and processed with the bundled Pix4DMapper software? Considering that UAV surveys for precision agriculture typically are multitemporal, the original question can be rephrased: what is the consistency between repeated surveys? In other words, when two datasets highlight differences for a certain part of a field, to what extent is this due to acquisition and processing errors, and to what extent does this point out a variation in the status of the crop? The importance of such questions is confirmed by the fact that only a few papers in the literature have explored them to date.

The present work studies geometric and radiometric consistency of two overlapping datasets, acquired with a Sequoia camera and processed with the bundled software. We focus on geometry to avoid the influence of its inconsistencies on the quality of the radiometry. A distinctive feature of the paper is that the geometric consistency is not assessed by means of a (generally limited) number of check points (CPs), as is usually done. Instead, we assess it by exhaustively evaluating the distance between the whole generated point clouds. We investigate radiometry as well, because it is the main source for agronomic studies. Moreover, we compare datasets acquired almost at the same time. This is a strength, as the difference assessed in vegetation indices can only be attributable to sensor noise, and possibly to issues in the radiometric calibration procedure.

#### **2. Methods**

#### *2.1. The Equipment*

The dataset was acquired with the HEXA-PRO™ UAV, which is operated by the Laboratory of Geomatics of the University of Pavia and is shown in Figure 1a. The vehicle was made by a small Italian company named Restart® and has the following main characteristics: 6 engines (290 W each one), Arducopter-compliant flight controller, maximum payload of 1.5 kg (partly used by the gimbal, weighting 0.3 kg), flight autonomy of approximately 15 min. The UAV was equipped with a Parrot Sequoia camera (see Figure 1c). Sequoia has a high-resolution RGB camera with a 4608 × 3456 pixel sensor, a pixel size of 1.34 μm, and a focal length of 4.88 mm; the ground sampling distance (GSD) is 1.9 cm at 70 m height above ground level (AGL). Sequoia also has four monochrome cameras that are sensitive to the following spectral bands: green (G, 530–570 nm), red (R, 640–680 nm), red-edge (RE, 730–740 nm), and near-infrared (NIR, 77–810 nm). Their resolution is 1280 × 960, with a pixel size of 3.75 μm and a focal length equal to 3.98 mm; the GSD is 6.8 cm at the 70 m flying height (AGL), which was adopted for the described survey.

**Figure 1.** The equipment operated by the Laboratory of Geomatics of the University of Pavia: (**a**) the HEXA-PRO™ unmanned aerial vehicle (UAV) used for the survey; (**b**) the Airinov calibration target supplied with the camera; (**c**) the Parrot Sequoia camera (the imaging and irradiance sensors are shown); (**d**) an example of the artificial markers used.

#### *2.2. The Block Structure*

On September 13, 2017, a photogrammetric survey was performed on the Santa Sofia farmstead, near Pavia, Northern Italy (Figure 2a). The test-site is a flat area totaling about 36 ha, used exclusively to cultivate rice. The whole acquisition was obtained by five flight missions, the planning for which is shown in Figure 2b, where the optical orthomosaic, which was used as background, was derived from a previous survey. In total, the project constituted about 1300 multispectral images, each composed of four bands. The AGL height was 70 m and image overlapping was 80% and 60% along- and across-track, respectively. The Sequoia camera was adopted, as previously mentioned. Twelve markers were placed on the ground and surveyed with the network real-time kinematic (NRTK) GPS mode. Virtual reference station (VRS) differential corrections were applied via connecting networked transport of RTCM (Radio Technical Commission for Maritime) via internet protocol (NTRIP - Networked Transport of RTCM via Internet Protocol) to the GNSS positioning service of "Regione Piemonte and Regione Lombardia" [34]. GCP coordinates have a planimetric and altimetric accuracy of 2–3 cm and 4–5 cm, respectively. GCPs were constituted by artificial markers with black and gray diamond shapes (Figure 1d); marker positions are illustrated in Figure 2b. At the beginning of each flight, the recommended radiometric calibration procedure was performed by acquiring the calibration target (Figure 1b).

The present paper will only focus on flights 3 and 4, as these had a methodological purpose. The overlapping area allowed us to deeply analyze geometric and radiometric congruency under several processing scenarios (as described in Section 2.3), because it is quite wide (23 ha) and encompasses 4 GCPs.

**Figure 2.** Unmanned aerial vehicle (UAV) survey framework: (**a**) site location; (**b**) the sub-block compositions. Note: light blue lines represent the flight outlines where the overlapping areas are clearly visible. The one considered in the paper is highlighted in red and includes four ground control points (GCPs), named 6, 7, 8, and 9. GCP locations are reported with green triangles. Coordinate reference system (CRS): WGS84/UTM 32N. Central coordinates (E, N): 506500, 5005600.

#### *2.3. The Photogrammetric Processing*

The photogrammetric project was carried out with Pix4Dmapper Pro, version 4.4.9. Only the four multispectral channels were considered, having 6.8 cm GSD; higher resolution RGB imagery was disregarded, as it is recorded in the JPEG format with a high compression factor, and has low quality compared to photogrammetry requirements. The processing followed the usual pipeline [35,36]: image alignment, tie point extraction, manual measurement of GCPs and CPs, bundle block adjustment (BBA), generation of dense point clouds, digital surface modeling (DSM), and creation of orthomosaic and reflectance maps. The software allows only one set of calibration target images to be used per project, so the photogrammetric processing followed a single-block approach. Four scenarios were depicted based on georeferentiation methodology and radiometric processing:


blocks with a large number of images and an overlapping area. It should ensure that radiometric differences caused by a misalignment in the dense point clouds are avoided. Scenario Ig/Jr was used only in radiometric assessment.

4. Joint georeferentiation/independent radiometric processing (Jg/Ir) scenario: the two blocks were jointly orientated, and the obtained exterior orientation parameters were then transferred to a single-block project for generation of dense point clouds and reflectance maps. In this scenario, possible radiometric inconsistencies due to separate computation of interior and exterior orientation parameters are eliminated. This scenario was used in both geometric and radiometric assessment.

The bundle block adjustment parameters were set according to the described scenario, since they differ in terms of calibration method and camera optimization. In DG scenario, the calibration method was set to the "alternative" option. This choice is recommended when the surveyed area is flat (as in this case) and there is availability of good image geolocation; for the Sequoia sensor, the used geolocation comes from the on-board GPS receiver, even if its quality is low, as discussed before. For camera optimization, external parameters were all re-estimated, while for the internal ones they were adopted from the camera model that is delivered by Sequoia directly into the EXIF (Exchangeable image file) section of each image. As we knew from the Pix4D technical support, the parameters delivered into the EXIF are individually determined for each item at the factory. Their reliability is good, as reported in [11], in which the changes between nominal and optimized camera parameters were as low as 0.01%. In Ig/Ir and Ig/Jr scenarios, the calibration method was again set to "alternative". For camera optimization, since the GCPs were imported and measured on each of the two blocks, both external and internal parameters were optimized. Finally, Jg/Ir is a two-step scenario in which the two blocks were jointly processed, and so the obtained internal and external parameters were used to separately generate the dense point clouds for each block. For the first step (image orientation), the parameters were set as equal to Ig/Ir; for the second step (single-block dense point cloud generation), the calibration method was set to geolocation-based, since accurate positioning and orientation are available from the first step. Besides, in this case, neither interior nor external parameters were optimized because they were directly imported in the first step of the project.

All dense point cloud generation was performed by adopting the default options: half image size resolution images, point density was set to optimal, and a cloud point was accepted only if it was positively matched in at least three images. The average density was between 11 to 14 points per m3. In a preliminary test, the original image size resolution was also evaluated, but higher point density did not significantly improve the generation of orthophotos and reflectance maps; the requirements for precise agriculture are lower in comparison to other applications, such as 3D mapping, and the obtained resolution was considered satisfactory for the research aims.

Pix4Dmapper allows generation of orthophotos and reflectance maps during step 3 of the processing, together with the computation of the DSM. In this study, products were generated with GSD equal to 0.10 m and project settings were the same for all considered scenarios. Reflectance maps were generated by setting camera and sun irradiance correction in the radiometric processing and calibration panel. This allows one to apply corrections to the camera parameters stored in the image metadata (i.e., vignetting, dark current, ISO), as well as for the sun irradiance information acquired with the proper sensor (see Figure 1c). Images of the calibration target are required to perform corrections. Hence, during the survey, the prescribed radiometric calibration procedure of the Parrot Sequoia camera was performed and the suitable calibration target (see Figure 1b) was imaged several times, with different exposure times. Acquisitions were taken at the beginning of each flight, so that different calibration data were stored for each flight, ensuring similar sky and illumination conditions between calibration images and flight images.

For the radiometric processing and calibration, calibration images with the highest value of exposure time were retained and the software automatically detected target on them, defining the proper reflectance values for each spectral band as equal to 0.172, 0.215, 0.266 and 0.369 for green, red, red-edge, and NIR, respectively.

#### *2.4. Geometric Consistency Assessment*

The iterative closest point (ICP) methodology was adopted to register overlapping point clouds, evaluate their distance, and estimate their geometric consistency. As is well-known in literature [37–39], ICP is a procedure aiming to align point clouds without requiring the identification of homologous points. It starts by associating each point of cloud A to its closest point belonging to cloud B, then a coordinate transformation (typically a roto-translation, having six parameters, also known as a rigid body transformation in literature [40]) is estimated, based on the obtained coupled points, and applied to one point set. The procedure is iterated until the latest estimated transformation is negligible. A dedicated Matlab procedure was specifically developed at the University of Pavia, implementing ICP and including some unique features. Procedure flowcharts are reported in Figure 3 and Figure 5.

**Figure 3.** Flowchart for data preparation.

Preliminarily (Figure 3), the common enveloped area for the two point clouds is determined. The clouds are then trimmed according to this area, adding a further precautionary buffer to avoid edge effects; the buffer value is set manually. Therefore, a subset of points is extracted (the so-called skeleton) from each of the point clouds that have a variable density. The skeleton is constituted by squared meshes with sides measuring 2 m and belonging to two classes. There are skeleton points inside meshes. Those lying on flat terrain contain 1 pt/m2, and therefore the spacing is 1 m. The others, which lie where there are ditches and escarpments, contain 64 pt/m2, with a spacing of 0.125 m. The skeleton was adopted to reduce the complexity of the calculation and to avoid flat terrain parts, where most of the original points are de-facto overweighed.

The classification of each mesh of the skeleton was performed by selecting all the original cloud points lying in the mesh and estimating the interpolating plane. By imposing suitable criteria concerning the residuals (low residuals mean flat areas) and the deviation from the vertical of the plane normal, the two classes (flat and variable terrain) were decided quite effectively. In the current version, the thresholds used for skeleton classification must be tuned by the operator and inserted manually; an improved one is under development, based on machine-learning. An example of the skeleton structure is reported in Figure 4 for sub-block 3, showing the skeleton points for flat and steep parts of the terrain.

**Figure 4.** The skeleton structure for sub-block 3. (**a**) Overview: the parts shown in light grey have a lower resolution of 1 pt/m2; the darker parts have a higher resolution of 64 pt/m2. (**b**) Detailed view: individual skeleton points are visible, and their different densities can be easily detected.

After defining the skeleton as described above, we used it to define a new point cloud. By this, we mean a set of points for which we know the 3D coordinates, the normal vector of the surface at their position, and the color, totaling 9 descriptors. Original point clouds A and B will be referred with the acronyms PCA and PCB, respectively, while point clouds obtained from skeletons will be named skeleton point clouds A and B, shortened to SPCA and SPCB, respectively. Each element of the skeleton point clouds has 9 descriptors, as stated before; some were known from the definition and some were calculated. For each skeleton mesh, the fitting plane was used to estimate the height of the interior skeleton points. The plane's cartesian equation was also used to obtain its normal vector. Finally, the color was determined for each skeleton point by picking that of the closest original point.

To perform quality assessment, data filtering, and further analysis, data used for each plane estimation was stored in a complex data structure. It is named cell array in the Matlab environment and can be thought of as a matrix where each element can store any kind of data structure. We created a cell array containing as many rows as the meshes constituting the skeleton and with 9 columns (the same number as the descriptors associated with each point cloud, but purely by chance—there is no relationship). For each mesh, the associated 9 cells contain the planimetric coordinates of the central point; mesh size; mesh point density, where we counted the number of the cloud points lying in the mesh and divide this by the area; mesh planimetric bounding box; mesh edge, intended as the points defining the corresponding polygon; coordinates of the original cloud points lying in the mesh; interpolating plane's normal vector; parameters of the plane and its fitting goodness (the so-called Matlab gof (goodness of fit) data structure); and the descriptors of the skeleton points located inside.

Once SPCA and SPCB were created, they were used to fit the ICP transformation. The procedure (Figure 5) is iterative; it starts with the original SPCA and SPCB and finally produces the parameters of a six-parameter 3D rigid transformation aligning SPCB to SPCA. Each iteration:


*Appl. Sci.* **2019**, *9*, 5314


**Figure 5.** Flowchart for ICP procedure. Note: K-D Tree, K-dimensional tree algorithm; MDL, KD tree model object; CTP, current transformation Parameters; TP, transformation parameters.

The process is stopped when the latest estimated transformation is negligible. The above-described process was coded in a Matlab toolbox coded at the University of Pavia. Among other features, the NLLS problem is solved with a robust approach, based on the Huber method. Moreover, the procedure takes advantage of the k-d tree functionalities (k-dimensional tree) to speed up point coupling, which are available in the used Matlab environment [40]. The k-d tree engine is trained for SPCA, which was kept fixed throughout the procedure.

The mathematical formulation [41,42] of the estimation for the coordinate transformation is

$$\text{CPTopt} = \arg\min\_{T} \sum \left( ( (\text{CPT} \cdot \text{SPC} B\_{i} - \text{SPCA}\_{i}) \cdot n\_{i} )^{2} \tag{1}$$

where *SPCBi* is the generic point of the skeleton B; *SPCAi* is the correspondent point of the skeleton A, derived by nearest neighbor searching; *ni* is the normal vector at point *SPCAi*; *CPT* is the 4 × 4 3D rigid-body transformation matrix estimated from previous iterations; and *CPTopt* is the 4 × 4 3D rigid-body transformation matrix estimated during the current iteration.

#### *2.5. Radiometric Consistency Assessment*

Radiometric consistency was assessed by computing, pixel by pixel, differences for the co-registered reflectance maps in the overlapping area of photogrammetric blocks 3 and 4. Respective statistics were also analyzed. Considering that Sequoia is a sensor mainly dedicated to agricultural applications, assessment was conducted also for some vegetation index (VI) maps, since they commonly represent a proxy of vegetation parameters to be used for agronomy purposes. VI maps were computed in Matlab, by applying an index formula to proper reflectance maps (Table 1).


**Table 1.** Vegetation indices (VIs) used in this study.

For Ig/Ir and Jg/Ir scenarios, maps derived from blocks 3 and 4 were directly compared, while Ig/Jr scenario was checked with respect to the single blocks of Ig/Ir scenario (see Section 2.3 for more details about scenario characteristics). From here on, maps are identified with the names "3 Ig/Ir", "4 Ig/Ir", "3 Jg/Ir", "4 Jg/Ir", "Ig/Jr", where "3" stands for block 3, and "4" for block 4. DG scenario was not considered for radiometric assessment.

Moreover, since no ground truth was available, the reliability of reflectance and VI maps was evaluated by comparing maps with the one obtained from Sentinel-2 (S2) imagery. Indeed, a Sentinel-2 acquisition two days after the survey (September 19, 2017) was available. Maps derived from the photogrammetric blocks (having a GSD equal to 0.10 m) were upscaled with a nearest-neighbor resampling to 10 m spatial resolution, to match Sentinel-2 imagery resolution. Correlation analysis was applied and statistics were performed on differences in terms of single bands and radiometric indices.

Although a comparison with ground truths calculated with a spectroradiometer would have been more effective, a test on compatibility between Sequoia and S2 data is also of scientific relevance, given the growing interest in the integration of data acquired from satellite and UAV platforms [48] for environmental applications [49,50], including PA [51–53], both from research and applied points of view.

#### **3. Results**

#### *3.1. Geometric Consistency*

#### 3.1.1. Reliability of the ICP-Derived Transformations

Since geometric consistency is based on ICP, it is mandatory to find a way to quantify the quality of this procedure, because if ICP fails, the estimated distance will not correspond to the actual one. It is known that the ICP estimation is not always reliable, especially when it is used to register almost flat clouds, as in our case. Registration is performed between two point clouds, namely PCA and PCB; it is possible to estimate the transformation PCA-to-PCB, which when applied to PCA, aligns it to PCB. This transformation is constituted by a roto-translation (i.e., the composition of a 3D shift

and a 3D rotation). Of course, it is possible to estimate the PCB-to-PCA transformation, which should coincide, aside from uncertainties, with the inverse of the first one. We used the comparison between the estimated and calculated inverse transformation to infer the precision of our estimations.

Results are shown in Table 2, where columns 3–5 report the components of the estimated shift (delta E, delta N, delta H) in meters and columns 6–8 show the rotation angles (ω, ϕ, κ) in degrees. Rows are grouped in fours, with each chunk being associated with one processing configuration. Row 1 reports the transformation (3-to-4) aligning point cloud 3 to cloud number 4. Row 2 shows the parameters of the calculated inverse transformation. Row 3 displays the parameters of the 4-to-3 estimated transformation and row 4 shows the differences between rows 2 and 3.


**Table 2.** Reliability of the iterative closest point (ICP)-estimated transformations to be used for cloud registration.

Scenario abbreviations: direct georeferencing (DG); independent georeferentiation/independent radiometric processing (Ig/Ir); joint georeferentiation/independent radiometric processing (Jg/Ir).

Excluding the DG scenario, in which direct georeferencing is adopted and point clouds are slightly deformed, the worst residual is 6 cm for the shift components and 0.06 deg for rotations. A distance-equivalent error (*e*) can be computed for the resulting angular residual (α) by assuming a distance (*d*) of 100 m, corresponding to the half-width of the considered test area. By applying the simple formula *e* = *d*α, where the angle α is expressed in radiants, it can be found that *e* = 12 cm.

Now, we must consider the granularity of the datasets (i.e., the points' linear spacing). For the skeleton, this is considered for ICP estimation; as already explained, the spacing is 12.5 cm for dense parts and 100 cm elsewhere. As residuals of the transformation equal the discretization size of the considered datasets, we consider the estimated transformations reliable and precise.

#### 3.1.2. Assessment of the Distance between Overlapping Blocks

There are three processing scenarios, and for each of them the distance between the two overlapping clouds (blocks 3 and 4) was assessed. Given two clouds, the ICP procedure was used to estimate the rigid transformation to register point cloud-A to point cloud-B. For the ICP procedure, the skeleton structures were used, as explained in Section 2.4, while for distance evaluation the original point clouds were considered. Once the transformation was applied, each point of B was coupled with the closest point from A; the components and the norm of the connecting vector were stored, together with the two indices addressing the selected point in the lists representing the two point-clouds. To work out the distance before ICP, the original point clouds were used with the same point couplings mentioned before.

Limited data cleaning was performed. First, a buffer (0.9 scale factor) was created along the border of the analyzed regions and points inside were disregarded. The goal was to ignore border effects, where the geometry of photogrammetric measurements is weak, and consequently model deformations might occur. Furthermore, a limited blunder rejection was performed. The empirical cumulative distribution function (CDF) of the 3D distances between coupled points was calculated and the pairs corresponding to values exceeding the 99th percentile were discarded.

Point cloud distances were assessed for the three scenarios, and descriptive statistics were applied to 4-tuples constituted by x, y, and z components of the displacement vector plus its norm. All the results are shown in Table 3. The first column reports the identifier of the scenario and how many point couples were used to evaluate the surface distance; columns 3 to 5 focus on the three components of the original clouds, while column 6 focuses on the 3D distance.


**Table 3.** Summary statistics of the 3D distance between overlapping point clouds.

The DG scenario shows large values, as expected, as the overlapping point clouds have an average 3D distance of 3.24 m. The considered scenario is based on direct georeferencing and the reported results confirm that the Sequoia's on-board GPS receiver is unfit for georeferencing photogrammetric products. This is not a surprise for the authors, nor should it be for any aware user. However, in times of widespread use of photogrammetry [54], we thought it was worth noting. Ig/Ir and Jg/Ir scenarios both adopt GCPs within different adjustment strategies. The RMSE values of the residuals for the single components *x*, *y*, and *z* are range between 6 and 11 cm. Considering the already mentioned granularity of the analyzed datasets, the reported figures highlight that the overlapping point clouds are optimally aligned and consistent.

Maps of the 3D distances between the overlapping point clouds are meaningful. Figure 6 shows them for all the three scenarios assessed. Remarkably, the three sub-figures shown adopt the same color map, even if Figure 6b,c only shows a small part of it. We also remark that values shown in Table 3 and plots reported in Figure 6 are related to the original point clouds, as they were generated by the photogrammetric procedure. ICP was only used to properly couple points belonging to different clouds in order to conveniently evaluate their distance.

Figure 6a highlights that in scenario DG, the clouds are quite far. Moreover, the map of the distances is a sort of a ramp, meaning that the two clouds are not simply displaced but are also affected by a significant rotation. The other two sub-figures are related to Ig/Ir and Jg/Ir scenarios and confirm that the distances are limited in size and are substantially constant. Moreover, Jg/Ir scenario shows lower values above the fields, where the terrain is flat, while distances are slightly greater beside dirt roads and ditches, where low vegetation is present.

**Figure 6.** The cloud 3D distance maps between areas 3 and 4, expressed in meters: (**a**) direct georeferencing (DG); (**b**) independent georeferentiation/independent radiometric processing (Ig/Ir); (**c**) joint georeferentiation/independent radiometric processing (Jg/Ir).

#### *3.2. Radiometric Consistency*

#### 3.2.1. Assessment of the Radiometric Differences between Overlapping Blocks

For the three processing scenarios, differences were calculated pixel by pixel among corresponding reflectance and VI maps in the overlapping area. While Ig/Ir and Jg/Ir scenarios were independently evaluated, Ig/Jr scenario was compared to the single blocks of Ig/Ir scenario (see Section 2.5). Descriptive statistics for differences calculated on reflectance maps are shown in Table 4, and results of VIs maps are reported in Table 5. Although differences have similar ranges, it is important to remember that reflectance maps have values in the range [0, 1], while values for VIs maps are in the range [−1, 1].

The computed RMSE values are quite close to zero for all cases, but significant differences among single reflectance maps and VI maps can be stressed, considering minimum and maximum absolute values. In particular, differences with maximum and minimum values above 0.4 are calculated for the NIR maps, differences reach values close to 0.3 for the red-edge map, and lower values are registered for the green and red maps, with minimum and maximum absolute values below 0.2 for the red maps in some cases. A similar behavior is also evident for the VI maps, where the differences calculated on NDVI maps assume lower RMSE values, while maximum and minimum values even greater than 0.5 are calculated for many VIs. The comparison between the statistics computed for Ig/Ir and Jg/Ir scenarios shows that both reflectance and VI map differences reach very similar values.


**Table 4.** Summary statistics of the differences between reflectance maps in the overlapping area.

**Table 5.** Summary statistics of the differences between VI maps in the overlapping area.


To assess the significance of the calculated values, the differences are presented in the form of box and whisker plots. Figure 7 reports box and whisker plots for differences computed on reflectance maps, while in Figure 8 VI results are shown. The plots do not refer to the Jg/Ir scenario, as similar results are obtained with respect to Ig/Ir scenario.

**Figure 7.** Box and whisker plots of differences computed on different reflectance maps in the overlapping area.

**Figure 8.** Box and whisker plots of differences computed on different VI maps in the overlapping area.

From the plots it is evident that results vary from map to map, but few general considerations can be drawn. Median values are overall around 0, while maximum and minimum values are outside of the confidence intervals and can be considered as outliers. For most cases, the variability of the differences is contained in the range [−0.2, 0.2]; thus, this interval of values is retained as significant for further analysis. The VIs can mitigate the effects of single reflectance maps, specifically the high differences registered for NIR maps are rather compensated in the NDVI maps. Moreover, with respect to the differences computed between single blocks (i.e., 3 Ig/Ir–4 Ig/Ir), results obtained considering Ig/Jr scenario (i.e., 3 Ig/Ir–Ig/Jr and 4 Ig/Ir–Ig/Jr) have narrower confidence intervals.

Spatial distribution of the differences in the overlapping area is shown in Figures 9–11. For the sake of brevity, only the most significative results are presented. As a matter of fact, similar results were registered for Ig/Ir and Jg/Ir scenarios. As regarding Ig/Jr scenario, green, NIR, and NDVI maps are shown, since the other maps have a similar spatial behavior. The remaining results are reported in Supplementary Materials (Figures S1–S3).

*Appl. Sci.* **2019**, *9*, 5314

**Figure 9.** Spatial distribution of differences in the overlapping area. Ig/Ir scenario: green (**a**), red (**b**), red-edge (**c**), NIR (**d**), NDVI (**e**), GNDVI (**f**), NDRE (**g**), NDVIre (**h**), and NGRDI (**i**).

**Figure 10.** Spatial distribution of differences in the overlapping area for Ig/Jr scenario, with respect to block 3: green (**a**), NIR (**b**), NDVI (**c**).

**Figure 11.** Spatial distribution of differences in the overlapping area for Ig/Jr scenario, with respect to block 4: green (**a**), NIR (**b**), NDVI (**c**).

A clear spatial pattern can be noted from the plots—the reflectance values tend to be overestimated as moving away from the center of the block (i.e., approaching the borders of the block); an analogous effect is visible in VI difference maps. This effect is more evident considering the differences calculated between the single blocks of Ig/Ir scenario (Figure 9). Tt is less evident when introducing also Ig/Jr scenario (Figures 10 and 11). No difference or very small differences are found in NDVI maps for all considered cases, which are uniformly distributed with no specific spatial profile in the overlapping area of blocks 3 and 4.

#### 3.2.2. Comparison with Sentinel-2 Imagery

As described in Section 2.5, the reliability of Sequoia maps was assessed with respect to Sentinel-2 data to evaluate the feasibility of data integration. First, an upscaling of maps derived from Sequoia imagery was required, then correlation analysis was computed (N = 265 samples). Results for the correlation analysis are reported in Figure 12 and map statistics are summarized in Table 6. For the sake of brevity, only results for NDVI are shown. As a matter of fact, other studies are present in the literature focusing on the comparison of NDVI only [51,52]. RMSE values reported in Table 6 were calculated by considering the NDVI map from S2 imagery as a reference.

**Figure 12.** Scatter plot and regression line for NDVI maps computed on S2 imagery with respect to Sequoia imagery: 3 Ig/Ir (**a**), 4 Ig/Ir (**b**), Ig/Jr (**c**). For each graph, the coefficient of determination (R2) and the Pearson's correlation coefficients () are reported (*p*-value <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16).

**Table 6.** Summary statistics of the NDVI maps computed from S2 imagery and Sequoia imagery, in the overlapping area.


The correlation with NDVI map from S2 imagery shows a good correspondence: coefficients of determination are 0.5197 for 3 Ig/Ir, 0.5249 for 4 Ig/Ir, and 0.4840 for Ig/Jr. The NDVI map with the highest correspondence against S2 imagery is the one derived from 4 Ig/Ir data, with Pearson's correlation coefficient and RMSE equal to 0.7245 and 0.0984, respectively. Nevertheless, the regression lines show a slight overestimate of Sequoia data compared to S2; NDVI maps from Sequoia imagery report higher values with respect to the S2 map (as also summarized by higher values for max. and mean in Table 6) and cover wider ranges (lower values for min. and higher values for STD in Table 6).

#### **4. Discussion**

#### *4.1. Geometric Consistency*

Geometric consistency was evaluated by computing the distance between the whole dense point clouds for the DG, Ig/Ir, and Jg/Ir scenarios. The Ig/Jr scenario was disregarded because although it is interesting for the study of radiometry, it is based on the extraction of a unique, common point cloud, and therefore there is nothing to check. Thanks to the high numerosity of the samples (the used clouds have between 10<sup>5</sup> and 10<sup>6</sup> points), we obtained very significant statistical figures. In addition, we could perform a sort of continuous evaluation of the distance between the overlapping surfaces.

The datasets considered are quite challenging for ICP, as the terrain is quite flat and the only height variations are due to medium-sized streams. Nevertheless, our developed algorithm proved to be reliable; indeed, we performed the closure check by composing direct and inverse transformations and obtained residual parameters (shift components and rotation angles), which are limited in size

and equivalent to the granularity of the considered point clouds, having an average distance of around 12 cm (see Section 3.1.1).

Cross checks of overlapping point clouds can be usefully applied when no or few CPs are available. Moreover, this allows detection of fine-grained deformations. Scenario DG is related to direct georeferencing, performed through the measurements of the Sequoia's on-board GPS receiver. Results are poor, not surprisingly, as the average distance between the clouds is above 3 m. The point sets are significantly shifted and rotated. This result was largely expected, but we think it is worth mentioning to warn newer photogrammetrists. Such georeferencing precision is not acceptable, even for precision agriculture applications. In Ig/Ir and Jg/Ir scenarios, bundle block adjustment was performed by means of 4 GCPs. Similar results are shown for both scenarios. Indeed, RMSE values of the components of the point-wise shift vectors are within 10 cm. It is the order of magnitude of the granularity of point clouds and of the pixel size of the generated orthomosaics. Such figures assure us that no difference or negligible radiometric differences are induced by geometric consistency issues.

In summary, two main findings come from geometry assessment: we have been able to reliably estimate a rigid 3D transformation by robust ICPs between almost flat point clouds; and we have demonstrated that the geometric consistency is good, so that the inconsistencies shown by radiometry have a different origin.

#### *4.2. Radiometric Consistency*

As already stressed by many authors in the literature [28–30,55], radiometric corrections are necessary when using sensors mounted on UAV for PA, but the ease of use and diffusion is limited. The radiometric processing available in Pix4Dmapper software for the Sequoia camera provides most of the corrections, including vignetting, dark current, exposure time, and sunlight irradiance, but omits other possible causes of radiometric variations [23]. First, this research points out that radiometric inconsistencies due to differences in the acquisition geometry remain unsolved. Reflectance values of pixels at the borders of the blocks tend to be overestimated, as a consequence of the inclination of the point of view during the photogrammetric survey. From Figure 9, it is evident that differences are not uniformly distributed, but present a clear spatial pattern. Higher difference values (absolute values) are measured at the borders of the overlapping area, while the lowest values approach the center of the area. This demonstrates the presence of a high edge effect on the reflectance maps, which must be considered during flights planning. In practical use, it is advisable to plan UAV surveys covering an area wider than the one of interest. Enlarging the survey area should guarantee uniformity in the acquisition geometry even in the edges, otherwise characterized by non-negligible radiometric distortions.

Radiometric differences are not affected by different geometric processing of the blocks, as confirmed by the similar values of the differences computed for Ig/Ir and Jg/Ir scenarios for both reflectance and VI maps (Tables 4 and 5). As can be noted from RMSE values reported in the tables, differences calculated for some VI maps are lower than values obtained for reflectance maps, meaning that some indices can decrease inconsistencies of single reflectance bands [56]. Considering the Ig/Jr scenario, the differences are moderate with respect to Ig/Ir scenario for both reflectance and VI maps, with mean values overall close to zero. The edge effect is also still evident from the spatial distributions shown in Figures 10 and 11, however with lower values, as is evident in the box plots in Figures 7 and 8. As a matter of fact, it should be recalled that the Ig/Jr scenario corresponds to the procedure recommended by Pix4Dmapper software to process large photogrammetric blocks. For adjacent blocks acquired with separate but temporally close flights, the recommended merging option can partially correct the effect of illumination geometry and mitigate radiometric inconsistencies in the overlapping areas between blocks.

There are still uncertainties regarding the obtained absolute values of reflectance and for the derived indices [56], and consequently in the quantitative use of the Sequoia data for the possible calculation of biophysical parameters. From the results reported in this study, it should be noted that in some areas, differences have values close or even larger than 0.2 (absolute value). Therefore, the different processing scenarios have an impact on the results in terms of radiometry. A difference of this magnitude cannot be neglected in the operational phase for precision agriculture applications; even more so if used for multitemporal surveys. As a matter of fact, the map that shows the most homogeneous values in all cases is NDVI, which is widely used in most agriculture applications [57–60].

Regarding the comparison with S2, which is limited in this paper to NDVI, it should be mentioned that despite the analysis being affected by the different geometric resolutions of sensors and acquisition platforms, a significative correlation is found between Sequoia and S2 maps. Following the approach of [52], the Pearson's correlation coefficient can be adopted as a map similarity measure. The obtained coefficients, which are close to 0.7, prove a coherence between the data collected from the different platforms and show similar spatial variability values of NDVI maps, which are to be interpreted as the same behavior in terms of crop vigor [61]. Therefore, the compatibility and integration of NDVI maps obtained by Sequoia and Pix4D systems should be feasible along with the Sentinel-2 products.

#### **5. Conclusions**

Even though producers and developers have made great efforts to enhance them, radiometric corrections leave significant radiometric distortions in orthomosaics obtained by Sequoia and Pix4D systems, which can result in biased absolute values. This study shows that relevant differences are found depending on flight geometry and block processing choices, with differences that can reach 20% of pixel values for single reflectance bands or VIs, thus reducing the effective use in PA. Moreover, available radiometric corrections do not guarantee uniform accuracy and consistency of results, and this can cause difficulties in comparing surveys carried on out different lightening conditions. Careful planning of the survey, together with proper choices of image processing, can enhance the results. Very high image overlap yields uniformity over a single block, and edge distortions can be reduced by surveying a wider area that includes the study area.

Nevertheless, for large surveys that imply the acquisition and processing of separated sub-blocks, the merge option suggested by Pix4D is effective in reducing radiometric inconsistencies in adjacent areas. This fact, together with the high correlation found with S2 products, proves that Sequoia is suitable for agronomic purposes, but great attention must be paid to the planning of the survey and to the data processing.

Therefore, it is necessary to increase the awareness in the use of sensors and semi-automatic data processing to deeply understand the strengths and weaknesses of UAV usage for PA. In this study, the choice to process the dataset following the proposed scenarios (Section 2.3) instead of a standard workflow was driven by the apparently impossibility of attributing the corresponding calibration set of images to each block. The new Sequoia+ sensor should bypass this issue because no calibration target is needed; imagery processing exploits a new fully automatic calibration pipeline in Pix4D. The authors do not have experience with this new camera release, however, they consider the proposed processing method an interesting and simple way to assess the performance of new sensors.

Finally, the applicability of the proposed method can be extended. In this paper, geometric and radiometric consistency was evaluated comparing results obtained from two almost contemporaneous flights, processed following a single-block approach. The same method can be used to evaluate consistency between two or more blocks acquired days or month apart; in other words, the method can be used to assess time-series.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/9/24/5314/s1. Figure S1: Spatial distribution of differences in the overlapping area for scenario Ig/Jr with respect to block 3: red (**a**), red-edge (**b**), GNDVI (**c**), NDRE (**d**), NDVIre (**e**), NGRDI (**f**). Figure S2: Spatial distribution of differences in the overlapping area for scenario Ig/Jr with respect to block 4: red (**a**), red-edge (**b**), GNDVI (**c**), NDRE (**d**), NDVIre (**e**), NGRDI (**f**). Figure S3: Spatial distribution of differences in the overlapping area. Scenario Jg/Ir: green (**a**), red (**b**), red-edge (**c**), NIR (**d**), NDVI (**e**), GNDVI (**f**), NDRE (**g**), NDVIre (**h**), NGRDI (**i**).

**Author Contributions:** Conceptualization, M.F. and V.C.; formal analysis, M.F., G.R., G.S., and V.C.; methodology, G.S. and V.C.; project administration, M.F.; software, M.F. and G.R.; supervision, G.S. and V.C.; visualization, M.F. and G.R.; writing—original draft, M.F., G.R., G.S., and V.C.; writing—review and editing, M.F., G.R., G.S., and V.C. **Acknowledgments:** Paolo Marchese is a technician of the Laboratory of Geomatics of the University of Pavia. He is the pilot of the UAV, plans the flight missions, carries out maintenance, and performs sensor integration; he is gratefully acknowledged here. We want to thank the farmhouse owner too, Roberto Manfredi, for always kindly allowing us to enter his property and carry out our activities.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
