Mapping Reflectance Anisotropy of a Potato Canopy Using Aerial Images Acquired with an Unmanned Aerial Vehicle

Viewing and illumination geometry has a strong influence on optical measurements of natural surfaces due to their anisotropic reflectance properties. Typically, cameras on-board unmanned aerial vehicles (UAVs) are affected by this because of their relatively large field of view (FOV) and thus large range of viewing angles. In this study, we investigated the magnitude of reflectance anisotropy effects in the 500–900 nm range, captured by a frame camera mounted on a UAV during a standard mapping flight. After orthorectification and georeferencing of the images collected by the camera, we calculated the viewing geometry of all observations of each georeferenced ground pixel, forming a dataset with multi-angular observations. We performed UAV flights on two days during the summer of 2016 over an experimental potato field where different zones in the field received different nitrogen fertilization treatments. These fertilization levels caused variation in potato plant growth and thereby differences in structural properties such as leaf area index (LAI) and canopy cover. We fitted the Rahman–Pinty–Verstraete (RPV) model through the multi-angular observations of each ground pixel to quantify, interpret, and visualize the anisotropy patterns in our study area. The Θ parameter of the RPV model, which controls the proportion of forward and backward scattering, showed strong correlation with canopy cover, where in general an increase in canopy cover resulted in a reduction of backward scattering intensity, indicating that reflectance anisotropy contains information on canopy structure. In this paper, we demonstrated that anisotropy data can be extracted from measurements using a frame camera, collected during a typical UAV mapping flight. Future research will focus on how to use the anisotropy signal as a source of information for estimation of physical vegetation properties.


Introduction
In optical remote sensing, the sensors are mostly passive systems and their measured signal therefore depends on the position and orientation of the sensor and the position of the sun relative to the observed surface. Differences in the intensity of the measured signal of natural surfaces as a function of observation and viewing geometry are the result of their anisotropic reflectance behavior. In its idealized form, this is commonly referred to as the Bidirectional Reflectance cover in a potato field where different nitrogen fertilization treatments were applied to sections in the field to study the effect of canopy development on reflectance anisotropy.

Study Area
The study area was a field with consumer potatoes (Solanum tuberosum L., cultivar Fontane) of approximately 450 × 200 m, located south of the village Reusel (51 • 59 47.9 N, 5 • 9 34.5 E) on the Dutch-Belgian border, in the province of Noord-Brabant (The Netherlands). The planting date was 16 April 2016. In this field, an experiment was performed to evaluate the effect of split-level fertilization on potato yield. Different initial fertilization levels (before planting) and sensor-based additional sensor-based fertilization (during crop growth) were applied to zones in the field ( Figure 1 and Table 1). The sensor-based variable-rate fertilization was applied based on crop reflectance measurements performed by Fritzmeier ISARIA crop sensors (Fritzmeier Umwelttechnik, Großhelfendorf, Germany) mounted in front of the tractor providing information on nitrogen need of the crop on a weekly basis. The variable rates for nitrogen were based on measured indices from the crop sensor and decision rules derived from previous dedicated fertilization experiments between 2010 and 2013 [43], and applied using a GNSS guided fertilizer spreader. Eight 30 × 30 m plots ( Figure 1: A-H) were created in the different fertilization zones where weekly measurements of leaf area index (LAI) and leaf chlorophyll content were performed to monitor the status of the potato plants. Sensor-based nitrogen fertilization was applied to half of the plots on 28 June 2016, 15 July 2016, and 9 August 2016. In addition, on 7 July 2016, fertilization with potassium of 60 kg/ha was applied to all plots.
Remote Sens. 2017, 9,417 3 of 23 LAI and canopy cover in a potato field where different nitrogen fertilization treatments were applied to sections in the field to study the effect of canopy development on reflectance anisotropy.

Study Area
The study area was a field with consumer potatoes (Solanum tuberosum L., cultivar Fontane) of approximately 450 × 200 m, located south of the village Reusel (51°59′47.9′′N, 5°9′34.5′′E) on the Dutch-Belgian border, in the province of Noord-Brabant (The Netherlands). The planting date was 16 April 2016. In this field, an experiment was performed to evaluate the effect of split-level fertilization on potato yield. Different initial fertilization levels (before planting) and sensor-based additional sensor-based fertilization (during crop growth) were applied to zones in the field ( Figure  1 and Table 1). The sensor-based variable-rate fertilization was applied based on crop reflectance measurements performed by Fritzmeier ISARIA crop sensors (Fritzmeier Umwelttechnik, Großhelfendorf, Germany) mounted in front of the tractor providing information on nitrogen need of the crop on a weekly basis. The variable rates for nitrogen were based on measured indices from the crop sensor and decision rules derived from previous dedicated fertilization experiments between 2010 and 2013 [43], and applied using a GNSS guided fertilizer spreader. Eight 30 × 30 m plots ( Figure  1: A-H) were created in the different fertilization zones where weekly measurements of leaf area index (LAI) and leaf chlorophyll content were performed to monitor the status of the potato plants. Sensor-based nitrogen fertilization was applied to half of the plots on 28 June 2016, 15 July 2016, and 9 August 2016. In addition, on 7 July 2016, fertilization with potassium of 60 kg/ha was applied to all plots.  Table 2). The base map is an RGB composite (R = 674 nm, G = 559 nm, and B = 500 nm), produced with the data collected during the flight.  Table 2). The base map is an RGB composite (R = 674 nm, G = 559 nm, and B = 500 nm), produced with the data collected during the flight. The ground-based apparent LAI (hereafter referred to LAI) was estimated based on weekly measurements with a Cropscan Multispectral Radiometer that used the weighted difference vegetation index (WDVI) [44] for estimating LAI. Calibration curves between the WDVI and LAI were determined by performing both Cropscan measurements and LAI2000 Plant Canopy Analyser measurements in the field at the same farm for potatoes grown in the years 2010-2014 [45]. In 2016, limited measurements with the LAI2000 were available, but validation with these measurements showed estimating LAI with Cropscan data was accurate. The percentage of canopy cover was estimated by classification of the pixels in orthophotos of the study area at high spatial resolution (~8 cm ground sampling distance (GSD)) into a vegetation and non-vegetation class. For this, we calculated the normalized difference vegetation index (NDVI) (Equation (1)).
where R 658 nm , R 674 nm , and R 848 nm are the reflectance factors at the subscripted wavelengths. We chose these bands because they are closest to Red Band 4 (665 nm) and NIR Band 8 (842 nm) of the Multi-Spectral Imager (MSI) onboard ESA's Sentinel-2 satellite, that are used for NDVI calculation [46]. The average of Rikola bands at 658 and 674 nm was calculated to approximate the reflectance of MSI Band 4 at 665 nm. Pixels with an NDVI > 0.7 were classified as vegetation. We chose this threshold value since it provided good results in separation of the vegetation and non-vegetation classes based on visual inspection.

UAV Flights
Two UAV flights, with an Altura AT8 octocopter, were performed during the growing season in the summer of 2016 (Table 2); one flight at an early growing phase (9 June 2016) and one flight at a late growing phase after the potato canopy was fully developed (19 July 2016). Between flights on the two dates, additional sensor-based fertilization was given to half of the plots ( Table 1). The flights were performed under clear-sky conditions to assure the strongest anisotropy effects [47]. The UAV was programmed to fly at 120 m above ground level with a horizontal speed of 4 m/s from waypoint to waypoint, resulting in four flight-lines that were spaced 20 m apart ( Figure 1). Images were collected continuously at an interval of 2.4 s. Flying at the aforementioned height and speed with the camera that we used in our study (Section 2.3) resulted in forward and sideways overlap of approximately 88% and 75%, respectively, with a GSD of approximately 8 cm.

Spectral Measurements
The sensor used for our measurements was a frame-based hyperspectral camera (Rikola Ltd., Oulu, Finland), which has a horizontal and vertical FOV of 36.5 • and collects images of 1010 × 1010 pixels. The camera uses the Fabry-Perot interferometer (FPI) technology to capture a series of snapshot images on programmatically selectable spectral bands in the visible-near infrared (VIS-NIR) wavelength range. An FPI system uses an optical filter with two semi-transparent mirror films facing each other. The wavelength of the transmitted light depends on the distance between these two films (the air gap). Thus, by taking a sequence of images while varying the distance between the films, a complete spectrum can be captured for each pixel. As the images on different wavelengths are collected at different times, each wavelength covers a slightly different extent and has a slightly different orientation due to the movement of the UAV during the collection of a full datacube. The number of spectral bands that can be sampled and the center wavelengths of those can be selected according to the user's preference. For our measurements, we sampled 16 bands in the 500-900 nm range (Table 3). A full datacube of 16 bands in raw format is 32.0 MB. The memory card that we used in our configuration, a 32.0 GB SanDisk CompactFlash memory card, has space for potentially storing 999 images of 16 bands. Prior to take-off for each flight, a dark current (image with lens cap on) measurement and image of a reflectance reference panel (gray 50% Spectralon panel (LabSphere Inc., North Sutton, NH, USA)) were taken. The HyperspectralImager 2.0 software by Rikola was used to convert the raw images to radiances. The radiance images were then transformed into reflectance factor images using the empirical line method and the measurement of the Spectralon panel. As the images were taken in clear atmospheric conditions from a relatively low altitude, performing advanced atmospheric correction to the images was not found to be necessary.

Orthorectification and Measurement Geometry
The reflectance factor images were aligned and ortho-and georectified in PhotoScan Professional 1.2.6 (AgiSoft LLC, St. Petersburg, Russia). Each of the 16 spectral bands had to be processed individually due to the movement of the UAV during the collection of a full datacube [48]. Per measured spectral band, the reflectance factor images were aligned (accuracy setting: highest) using the GPS data collected by the Rikola sensor during the flights as a reference. Thereafter, all bands were merged, and 12 natural ground control points (GCPs) (Figure 1) were used for accurate georeferencing of the data. The GCPs were manually selected in each of the 16 spectral bands. Unfortunately, there were no dedicated Real Time Kinematic (RTK) or GPS/GNSS ground control markers available at the time of the flights. The imagery was used to calculate and export a Digital Surface Model (DSM) of the target area at a 5 m ground pixel size in WGS84 UTM 31N coordinate system. Each reflectance factor image was georectified and resampled separately to the same 5 m ground pixels and exported as separate GeoTIFFs. We specifically chose a 5 m ground pixel size, because a pixel of this size is large enough to capture a representative proportion of potato rows and exposed soil between the potato rows, while it is small enough to capture the spatial variation caused by the different fertilization treatments applied to the plots in the field. Finally, the photogrammetric positions of the camera (x cam , y cam , z cam ) during each image acquisition were exported as an ASCII file. The PhotoScan processing was performed on a computer with a Windows 7 Enterprise (64 bit) operating system, with an Intel(R) Xeon(R) CPU E5-1650 0 processor at 3.20 GHz (12 CPUs), an AMD FirePro V5900 2 GB graphics card, and 8.00 GB of installed memory. The processing time of the aforementioned steps was approximately 3.5 h per flight.
After the PhotoScan processing, the location of the ground pixels x px , y px , z px were known. Based on the pixel and camera locations, we calculated the View Azimuth Angle (VAA) and View Zenith Angle (VZA) for each ground pixel and reflectance factor image combination. The solar illumination angle for each pixel, relative to the normal of the surface (β x px , y px ) was calculated using cos β x px , y px = cos θ s cos θ n x px , y px + sin θ s sin θ n x px , y px cos ϕ s − ϕ n x px , y px (2) where θ s and ϕ s are the sun zenith and sun azimuth angle (Table 2), respectively. θ n and ϕ n are the slope and aspect of the pixel, respectively, which were calculated from the DSM. The VZA (θ v ) was calculated using The VAA (ϕ v ) of each pixel was calculated using Equation (3) assumes that the measured surface is flat. To correct for variations in slope due to height differences in the field, Equations (2)-(4) were substituted in Equation (2), resulting in cos θ v corrected x px , y px = cosθ v x px , y px cosθ n x px , y px + sinθ v x px , y px sinθ n x px , y px cos ϕ v x px , y px − ϕ n x px , y px .
Throughout the rest of this article, cos θ v corrected is referred to as VZA. By applying the Equations (2)-(5), we were able to obtain the observation and illumination geometry of each individual pixel of each image at each spectral band (Figure 2). Besides a spatial location and reflectance factor (Figure 2b), each pixel now has a VAA and VZA (Figure 2c,d). Pixels that were captured by multiple cameras thus have multiple observation angles as well.

Data Analysis and Visualization
After determining the observation geometry for all ground pixels and aerial image combinations, we fitted the RPV model [42] to each ground pixel in the study area. The RPV parameters obtained for each pixel were used to create parameterized anisotropy maps of the study area and to analyze the anisotropy patterns that were observed in the different experimental plots. We choose the RPV model because it has been successfully used in a wide variety of anisotropy studies at different spatial scales, ranging from a centimeter scale in laboratory studies [49,50] to a kilometer scale in space-borne studies [51]. The RPV model defines the reflectance anisotropy with four parameters where the parameter controls the amplitude, the parameter controls the bowl ( < 1)/bell ( > 1) shape of the anisotropy curve, the parameter controls forward/backward scattering, and the parameter controls the hotspot effect around the direct backscattering angle. Mathematically the RPV model is defined as where

Data Analysis and Visualization
After determining the observation geometry for all ground pixels and aerial image combinations, we fitted the RPV model [42] to each ground pixel in the study area. The RPV parameters obtained for each pixel were used to create parameterized anisotropy maps of the study area and to analyze the anisotropy patterns that were observed in the different experimental plots. We choose the RPV model because it has been successfully used in a wide variety of anisotropy studies at different spatial scales, ranging from a centimeter scale in laboratory studies [49,50] to a kilometer scale in space-borne studies [51]. The RPV model defines the reflectance anisotropy with four parameters where the ρ 0 parameter controls the amplitude, the k parameter controls the bowl (k < 1)/bell (k > 1) shape of the anisotropy curve, the Θ parameter controls forward/backward scattering, and the ρ c parameter controls the hotspot effect around the direct backscattering angle. Mathematically the RPV model is defined as Remote Sens. 2017, 9, 417 8 of 23 where To fit the model to each ground pixel dataset, we applied an iterative optimization method that uses the Levenberg-Marquardt algorithm to minimize the sum of squared residuals (i.e., the difference between the measured reflectance factor at all observation angles and the reflectance factors modeled by the RPV model at the same observation angles) to estimate the parameters. During the fitting procedure, the range of possible solutions for the Θ parameter was limited between −1 and 1, the ranges of the ρ 0 and k parameter were left free. The hotspot parameter ρ c was fixed at 1, which disables the hot spot parameter, since we acquired no observations close to the hotspot.

Crop Development
In the data collected on 9 June 2016, the different initial fertilization levels were clearly visible as horizontal zones in NDVI map of the potato field ( Figure 3a). Plot C and Plot D, the plots in the zone that did not receive any initial fertilization, stood out with clearly lower NDVI values. During the period between the two measurement days, there were two moments at which the potato plants on Plots B, D, F, and H received additional sensor-based nitrogen fertilization. This resulted in higher NDVI values in the plots in this section compared to the plots in the section that did not receive any additional fertilization. The effect of the different additional fertilization levels on the NDVI is most clearly visible, as the difference between Plots C and D on 19 July (Figure 3b), where Plot D has clearly higher NDVI values. Figure 3c shows the development of the LAI based on Cropscan measurements throughout the growing season. During Flight 1 on 9 June 2016, the potato crop was still at an early stage of growth with a haulm (leaves and stems) length of 0.4-0.6 m. The potato plants were standing upright above the ridges, showing a clear row structure. The LAI of all plots increased until the last week of June ( Figure 3c). During the first two weeks of July, the LAI decreased, which can be explained by growth-reducing effects caused by a dry period with high temperatures in the area where the study site was located. The insufficient water availability due to the lack of rain, combined with the high temperatures can affect the stomatal capacity of leaves, which in turn can cause a reduction biomass, and thus LAI [52]. Flight 2 on 19 July 2016 occurred after the potato plants were flowered, and the haulm length increased to 0.6-0.9 m. Due to the drought and high temperatures, stems had fallen over and started filling up the space between the potato rows. This resulted in an increase in canopy cover at the date of the second flight, while the LAI remained more or less equal, or for some plots became even lower (Table 4). After the second flight, the improved water availability due to rainfall and combined with the fertilizer addition of July 15 (Table 1) resulted in an increase in LAI with the largest long-term effect for the Plots B, D, F and H (Figure 3c).

View Angle Coverage
Based on Equations (2)-(5), we determined the VAAs and VZAs of each pixel that was georeferenced during the Agisoft PhotoScan processing. Pixels that were captured from multiple camera positions were thus also viewed from multiple observation geometries. As an example, Figure 4a shows from how many camera positions each pixel in the study area was captured at 658 nm during the flight on 9 June 2016. An example pixel (the white square in between Plot C and Plot D) was captured from 32 different camera positions in all four flight lines (Figure 4b).

View angle coverage
Based on Equations (2)-(5), we determined the VAAs and VZAs of each pixel that was georeferenced during the Agisoft PhotoScan processing. Pixels that were captured from multiple camera positions were thus also viewed from multiple observation geometries. As an example, Figure  4a shows from how many camera positions each pixel in the study area was captured at 658 nm during the flight on 9 June 2016. An example pixel (the white square in between Plot C and Plot D) was captured from 32 different camera positions in all four flight lines (Figure 4b).  The nadir normalized reflectance (anisotropy factors (ANIFs) [3] observed from the different camera positions at 658 nm for the example pixel indicated a strong anisotropy effect, where the highest ANIFs were observed in the backward scattering direction in the solar principal plane (Figure 5a). This backward scattering anisotropy was weaker for the same pixel at 848 nm (Figure 5b). Due to the movement of the UAV during collection of a full datacube, the camera positions were slightly different for both bands. In addition, the example pixel was not captured from all the same camera positions at both wavelengths. For example, the example pixel was not captured from Camera Positions 12 and 147 at 848 nm. On the contrary, the example pixel was not captured from Camera Position 58 at 658 nm, while it was at 848 nm (compare Figure 5a,b).
Remote Sens. 2017, 9, 417 11 of 23 highest ANIFs were observed in the backward scattering direction in the solar principal plane ( Figure  5a). This backward scattering anisotropy was weaker for the same pixel at 848 nm (Figure 5b). Due to the movement of the UAV during collection of a full datacube, the camera positions were slightly different for both bands. In addition, the example pixel was not captured from all the same camera positions at both wavelengths. For example, the example pixel was not captured from Camera Positions 12 and 147 at 848 nm. On the contrary, the example pixel was not captured from Camera Position 58 at 658 nm, while it was at 848 nm (compare Figure 5a,b).

Anisotropy Maps
For every georeferenced ground pixel in the study area at each measured wavelength, multiangular observations as displayed in Figure 5 were available. By fitting the RPV model through these multi-angular observations at every pixel, we created anisotropy maps that show the spatial distribution of the model parameters. Figure 6 shows the parameter maps of the , , and parameters at 658 nm and at 848 nm on 9 June 2016. The parameter, which is closely related to the nadir reflectance, follows a typical vegetation trend, with low values at 658 nm ( Figure 6a) and higher values at 848 nm (Figure 6d). Plots C and D, which did not receive any initial fertilization and therefore had a lower canopy cover, showed clearly lower values at 848 nm. In general, pixels in both bands had < 1, indicating a bowl-shaped anisotropy pattern (Figure 6c,e). Both wavelengths showed clear backward scattering anisotropy ( < 0) over the full study area. This backward scattering was particularly strong at 658 nm ( Figure 6c) and less pronounced at 848 nm (Figure 6f).

Anisotropy Maps
For every georeferenced ground pixel in the study area at each measured wavelength, multi-angular observations as displayed in Figure 5 were available. By fitting the RPV model through these multi-angular observations at every pixel, we created anisotropy maps that show the spatial distribution of the model parameters. Figure 6 shows the parameter maps of the ρ 0 , k, and Θ parameters at 658 nm and at 848 nm on 9 June 2016. The ρ 0 parameter, which is closely related to the nadir reflectance, follows a typical vegetation trend, with low values at 658 nm ( Figure 6a) and higher values at 848 nm (Figure 6d). Plots C and D, which did not receive any initial fertilization and therefore had a lower canopy cover, showed clearly lower ρ 0 values at 848 nm. In general, pixels in both bands had k < 1, indicating a bowl-shaped anisotropy pattern (Figure 6c,e). Both wavelengths showed clear backward scattering anisotropy (Θ < 0) over the full study area. This backward scattering was particularly strong at 658 nm ( Figure 6c) and less pronounced at 848 nm (Figure 6f). On 19 July 2016, the potato canopy was fully developed, resulting in general in higher values at 848 nm (Figure 7d). The pixels of Plot C, which was the plot in the zone that did not receive any initial fertilization nor additional sensor-based fertilization, had clearly lower values at 848 nm than the other plots. The parameter remained more or less the same as on the first measuring day. The parameter, on the contrary, showed a strong increase at both wavelengths (i.e., the values became less negative), indicating a reduction of backward scattering intensity (compare Figures 6c,f and 7c,f). On 19 July 2016, the potato canopy was fully developed, resulting in general in higher ρ 0 values at 848 nm (Figure 7d). The pixels of Plot C, which was the plot in the zone that did not receive any initial fertilization nor additional sensor-based fertilization, had clearly lower ρ 0 values at 848 nm than the other plots. The k parameter remained more or less the same as on the first measuring day. The Θ parameter, on the contrary, showed a strong increase at both wavelengths (i.e., the values became less negative), indicating a reduction of backward scattering intensity (compare Figure 6c Figure 8 shows the interpolated polar plots of the ANIFs at 658 nm and 848 nm for all pixels within Plot C and Plot E, two plots that varied strongly in LAI and canopy cover. Pixels within a 3 m distance of tractor tracks, which were running through the center of the plots (Figure 1), were excluded from this analysis. At both wavelengths, there was a clear backward scattering anisotropy in the plots, which was most pronounced at 658 nm. During the flight on 9 June 2016, we observed a weaker backward scattering anisotropy in Plot C at 658 nm compared to Plot E, as can be observed by the higher density of isolines. This was likely due to the fact that the potato plant rows in Plot C were not continuous due to gaps in the potato rows, either because of missing plants in the rows or because of the small size of the plants (Figure 4b). This in turn resulted in fewer strong shadows and thereby a weaker backward scattering intensity. The rows of potato plants in Plot E, on the other hand, had fewer gaps and therefore displayed stronger shadows, and thus displayed a stronger backward scattering intensity. During the flight on 19 July 2016, the canopy of Plot E had fully developed into a closed surface with a row structure that was at that point hardly visible. This  Figure 8 shows the interpolated polar plots of the ANIFs at 658 nm and 848 nm for all pixels within Plot C and Plot E, two plots that varied strongly in LAI and canopy cover. Pixels within a 3 m distance of tractor tracks, which were running through the center of the plots (Figure 1), were excluded from this analysis. At both wavelengths, there was a clear backward scattering anisotropy in the plots, which was most pronounced at 658 nm. During the flight on 9 June 2016, we observed a weaker backward scattering anisotropy in Plot C at 658 nm compared to Plot E, as can be observed by the higher density of isolines. This was likely due to the fact that the potato plant rows in Plot C were not continuous due to gaps in the potato rows, either because of missing plants in the rows or because of the small size of the plants (Figure 4b). This in turn resulted in fewer strong shadows and thereby a weaker backward scattering intensity. The rows of potato plants in Plot E, on the other hand, had fewer gaps and therefore displayed stronger shadows, and thus displayed a stronger backward scattering intensity. During the flight on 19 July 2016, the canopy of Plot E had fully developed into a closed surface with a row structure that was at that point hardly visible. This resulted in less pronounced shadows cast by the potato rows and thus in a weaker backward scattering on this date. On both dates, the backward scattering intensity for both plots at 848 nm remained more or less similar.

Plot Statistics
Remote Sens. 2017, 9,417 14 of 23 resulted in less pronounced shadows cast by the potato rows and thus in a weaker backward scattering on this date. On both dates, the backward scattering intensity for both plots at 848 nm remained more or less similar. The average RPV parameters, based on the individually obtained RPV parameter values for every pixel within the experimental plots, are shown in Figure 9. Again, pixels within a 3 m radius of the tractor tracks were excluded. The parameter followed the pattern of a typical vegetation spectrum. In general, a higher value was observed in the green and NIR wavelength region for plots that had a higher canopy cover on 19 July 2016. This is most clear for the plots where there was a strong increase in canopy cover between the two dates ( Table 4). The parameter was in general <1 at all wavelengths, indicating a bowl-shaped anisotropy pattern for the whole sampled wavelength region. The parameter, like the parameter, followed in shape a vegetation spectrum. The strongest backward scattering effects (the lowest values) were observed in the visible wavelength region, with a minimum at 674 nm, where shadow effects are strongest and multiple scattering effects are absent due to the absorption of radiance by chlorophyll in this region. In the NIR region, where shadow effects are diminished due to high reflectance and transmittance by vegetation, the values were in general less negative than in the visible region. The values obtained for the flight on Day 2 were consistently higher (less negative) than the values obtained for the flight on Day 1, indicating a reduced backward scatting intensity with increased canopy cover. The average RPV parameters, based on the individually obtained RPV parameter values for every pixel within the experimental plots, are shown in Figure 9. Again, pixels within a 3 m radius of the tractor tracks were excluded. The ρ 0 parameter followed the pattern of a typical vegetation spectrum. In general, a higher ρ 0 value was observed in the green and NIR wavelength region for plots that had a higher canopy cover on 19 July 2016. This is most clear for the plots where there was a strong increase in canopy cover between the two dates ( Table 4). The k parameter was in general <1 at all wavelengths, indicating a bowl-shaped anisotropy pattern for the whole sampled wavelength region. The Θ parameter, like the ρ 0 parameter, followed in shape a vegetation spectrum. The strongest backward scattering effects (the lowest Θ values) were observed in the visible wavelength region, with a minimum at 674 nm, where shadow effects are strongest and multiple scattering effects are absent due to the absorption of radiance by chlorophyll in this region. In the NIR region, where shadow effects are diminished due to high reflectance and transmittance by vegetation, the Θ values were in general less negative than in the visible region. The Θ values obtained for the flight on Day 2 were consistently higher (less negative) than the Θ values obtained for the flight on Day 1, indicating a reduced backward scatting intensity with increased canopy cover.  Figure 10 shows the relation between the average RPV parameters of all pixels within each experimental plot and canopy cover at 658 nm and 848 nm, respectively. On both dates, an increased canopy cover resulted in a decrease in the parameter at 658 nm and in an increase at 848 nm. This trend can be explained by the lower reflectance of vegetation at 658 nm compared to the soil background reflectance at this wavelength. A higher canopy cover therefore resulted in higher reflectance, and thus a higher value. The opposite holds at 848 nm, where the soil reflectance was lower than the vegetation reflectance. On 9 June 2016 at 658 nm and on 19 July 2016 at 848 nm, the relations between and canopy cover were strong, indicated by the R 2 of 0.764 and 0.753, respectively.

RPV Parameters vs. Canopy Cover and LAI
The relation between the parameter and canopy cover was less obvious. On Day 1 at 658 nm, we observed a slight decrease in the parameter with increased canopy cover (R 2 = 0.478). On both dates at 848 nm, the relation between the parameter and canopy cover was very weak, indicated by the R 2 of 0.262 and 0.211 on Day 1 and Day 2, respectively.
On Day 1, we observed a strong decrease in the parameter (which indicates increased backward scattering intensity) with an increase in the canopy cover at 658 nm. On Day 2, on the contrary, we observed an increase in the parameter with increasing canopy cover. The decrease in the parameter on Day 1 happened most likely due to the aforementioned gaps in the potato rows that were present in the plots with low canopy cover on this day: the gaps in the rows result in weaker shadows and thus weaker backward scattering. On Day 2, when there were hardly any gaps left in the rows, we observed a decrease in backward scattering intensity with increasing canopy cover. This is likely due to the disappearance of the row-structure since the potato rows started growing into  Figure 10 shows the relation between the average RPV parameters of all pixels within each experimental plot and canopy cover at 658 nm and 848 nm, respectively. On both dates, an increased canopy cover resulted in a decrease in the ρ 0 parameter at 658 nm and in an increase at 848 nm. This trend can be explained by the lower reflectance of vegetation at 658 nm compared to the soil background reflectance at this wavelength. A higher canopy cover therefore resulted in higher reflectance, and thus a higher ρ 0 value. The opposite holds at 848 nm, where the soil reflectance was lower than the vegetation reflectance. On 9 June 2016 at 658 nm and on 19 July 2016 at 848 nm, the relations between ρ 0 and canopy cover were strong, indicated by the R 2 of 0.764 and 0.753, respectively.

RPV Parameters vs. Canopy Cover and LAI
The relation between the k parameter and canopy cover was less obvious. On Day 1 at 658 nm, we observed a slight decrease in the k parameter with increased canopy cover (R 2 = 0.478). On both dates at 848 nm, the relation between the k parameter and canopy cover was very weak, indicated by the R 2 of 0.262 and 0.211 on Day 1 and Day 2, respectively.
On Day 1, we observed a strong decrease in the Θ parameter (which indicates increased backward scattering intensity) with an increase in the canopy cover at 658 nm. On Day 2, on the contrary, we observed an increase in the Θ parameter with increasing canopy cover. The decrease in the Θ parameter on Day 1 happened most likely due to the aforementioned gaps in the potato rows that were present in the plots with low canopy cover on this day: the gaps in the rows result in weaker shadows and thus weaker backward scattering. On Day 2, when there were hardly any gaps left in the rows, we observed a decrease in backward scattering intensity with increasing canopy cover. This is likely due to the disappearance of the row-structure since the potato rows started growing into each other, narrowing and disappearing the space between the rows, which cause strong shadow effects. each other, narrowing and disappearing the space between the rows, which cause strong shadow effects. Finally, we studied the correlation of the RPV parameters with the canopy cover and LAI, respectively, for the different experimental plots at all measured spectral bands ( Figure 11). We calculated the Kendall's tau ranking correlation coefficient, since not all canopy cover, LAI, and RPV parameter values were normally distributed, and the number of plots and thus number of observations (n = 8) was rather low. The Kendall's tau correlation coefficient takes values between −1 and +1, where a positive value indicates that the ranks of both variables are increasing and a negative value indicates that the ranks of both variables are decreasing. The closer Kendall's tau gets to +1 or −1, the stronger the correlation between the variables. Since the flight on Day 1 was performed at an early phase of the potato crop before it reached maximum LAI and the flight on Day 2 at a late phase after reaching maximum LAI, both flights are first analyzed separately and then combined.
The correlation of the parameter with canopy cover (Figure 11a, upper graph) on Day 1 was negative in the visible part of the spectrum, whereas it was positive in the NIR. On Day 2, this correlation was only negative in the red. A significant correlation at 5% confidence level (p = 0.05) Figure 10. The relation between the ρ 0 , k, and Θ parameters of the RPV model of all pixels in each experimental plot (indicated by letters A-H) and the canopy cover, at 658 nm (black circles) and 848 nm (gray triangles), respectively, for both dates. The error bars indicate the standard deviation.
Finally, we studied the correlation of the RPV parameters with the canopy cover and LAI, respectively, for the different experimental plots at all measured spectral bands ( Figure 11). We calculated the Kendall's tau ranking correlation coefficient, since not all canopy cover, LAI, and RPV parameter values were normally distributed, and the number of plots and thus number of observations (n = 8) was rather low. The Kendall's tau correlation coefficient takes values between −1 and +1, where a positive value indicates that the ranks of both variables are increasing and a negative value indicates that the ranks of both variables are decreasing. The closer Kendall's tau gets to +1 or −1, the stronger the correlation between the variables. Since the flight on Day 1 was performed at an early phase of the potato crop before it reached maximum LAI and the flight on Day 2 at a late phase after reaching maximum LAI, both flights are first analyzed separately and then combined. was only observed in the NIR for both dates. The pattern for the correlation between and LAI was quite similar to the one found for canopy cover on both dates (compare Figure 11b with Figure 11a, upper graphs). Figure 11. Ranking correlation (Kendall's tau) between the , , and parameters of the RPV model and the canopy cover (a) and Cropscan LAI (b) for the flights of Day 1 and Day 2 separately and for the flights of both days combined. The gray-shaded areas indicate the significance levels: values above or below the dark-gray areas were significant for the analysis of both dates (n = 16) and values above or below the light-gray areas were significant for the analysis of the separate days (n = 8), both at the 5% confidence level.
The correlation of the parameter with canopy cover (Figure 11a, middle graph) was strongest (negative) on Day 1 in the visible part of the spectrum. Only the correlation at 658 nm was significant. The correlations between the parameter and LAI (Figure 11b, middle graph) were quite similar to the one for canopy cover. No correlations were significant at the 5% level.
The correlation between the parameter and canopy cover was negative in the visible part of the spectrum and positive in the NIR on Day 1, but they were not significant (Figure 11a, lower graph). On Day 2, the correlations were close to zero, except for positive, non-significant correlations in the visible part. The correlations between the parameter and LAI were similar to those of Figure 11. Ranking correlation (Kendall's tau) between the ρ 0 , k, and Θ parameters of the RPV model and the canopy cover (a) and Cropscan LAI (b) for the flights of Day 1 and Day 2 separately and for the flights of both days combined. The gray-shaded areas indicate the significance levels: values above or below the dark-gray areas were significant for the analysis of both dates (n = 16) and values above or below the light-gray areas were significant for the analysis of the separate days (n = 8), both at the 5% confidence level.
The correlation of the ρ 0 parameter with canopy cover (Figure 11a, upper graph) on Day 1 was negative in the visible part of the spectrum, whereas it was positive in the NIR. On Day 2, this correlation was only negative in the red. A significant correlation at 5% confidence level (p = 0.05) was only observed in the NIR for both dates. The pattern for the correlation between ρ 0 and LAI was quite similar to the one found for canopy cover on both dates (compare Figure 11b with Figure 11a, upper graphs).
The correlation of the k parameter with canopy cover (Figure 11a, middle graph) was strongest (negative) on Day 1 in the visible part of the spectrum. Only the correlation at 658 nm was significant. The correlations between the k parameter and LAI (Figure 11b, middle graph) were quite similar to the one for canopy cover. No correlations were significant at the 5% level.
The correlation between the Θ parameter and canopy cover was negative in the visible part of the spectrum and positive in the NIR on Day 1, but they were not significant (Figure 11a, lower graph). On Day 2, the correlations were close to zero, except for positive, non-significant correlations in the visible part. The correlations between the Θ parameter and LAI were similar to those of canopy cover. For LAI, the correlations were significant in the visible region on Day 1 (Figure 11b, lower graph).
When we combine the data collected on both days, the correlation between the ρ 0 parameter and canopy cover increased slightly over the whole spectrum (Figure 11a, upper graph). On the contrary, the correlation between the ρ 0 parameter and LAI decreased when both days were combined, which indicates that the ρ 0 parameter was more sensitive to canopy cover than to LAI, and differences between the structure of the potato crop on both dates were not well represented by this parameter. This was likely due to the fact that the LAI did not change much between the two dates, while there was a strong increase in canopy cover (Table 4). Over the whole spectrum, there was a positive correlation between the Θ parameter and canopy cover when the data of both dates were combined (Figure 11a, lower graph), suggesting that in general there was a decrease in backward scattering intensity (increase in the Θ parameter) when the canopy cover increased.

Discussion
In this paper, we demonstrated how to extract multi-angular observations from overlapping pixels captured by a frame-based camera during a typical UAV mapping flight. After orthorectification and georeferencing of the individual images taken at all measured spectral bands, we were able to calculate the observation geometry of all ground referenced pixels based on the pixel location and camera positions from which the pixels were observed.
During the two flights that we performed during the summer of 2016 over an experimental potato field, pixels were captured in up to 40 different images collected from different positions of the flight path of the UAV (Figure 4), resulting in up to 40 different observation geometries per pixel. Although the angular sampling coverage of these observations was limited, strong anisotropy effects occurred (Figure 5a,b). This highlights the importance of taking reflectance anisotropy into account when data is collected at off-nadir viewing angles. Either off-nadir data needs to be corrected and/or normalized to a standardized observation geometry [53], or the exact illumination and observation geometry needs to be considered in further analysis of the data [54].
In this study, we sampled a relatively large area, which resulted in a relatively low angular sampling coverage. To achieve a higher sampling density, we suggested focusing on smaller areas and design flight paths with higher overlap and/or to introduce flight lines in multiple azimuth directions covering the area of interest, like for example done by Honkavaara et al. [55] (p. 78, Figure 4). Our focus, however, was to derive anisotropy information during a typical UAV mapping flight, and we therefore did not test any alternative flight patterns than the one displayed in Figure 1. Another option to sample larger VZAs would be to mount the camera tilted under the UAV. However, this could complicate photogrammetric processing steps.
Several other methods to perform multi-angular reflectance measurements of vegetation targets using UAVs have recently been proposed ranging from spectrometers controlled by a gimbal [34,35] to the exploitation of the wide field of view of a hyperspectral line scanner [36]. Using these methods, clear anisotropy effects were observed for different vegetation targets at off-nadir viewing angles, where typically higher reflectance factors were observed in the backward scattering direction, with emphasis on view angles close to the hotspot position [34,36]. The drawback of these methods is that the measurements only provide the anisotropy information of a small area or the average anisotropy of a larger surface, respectively. Both of these methods lack information on the spatial distribution of the observed anisotropy effects. The method used in this paper does provide spatially explicit information on anisotropy. By fitting BRDF models, like we did with the RPV model, information on anisotropy can be stored and visualized in maps.
In this study, we used an FPI frame-based camera, which complicates the photogrammetric processing due to the movement of the UAV during the collection of a datacube [48]. In addition, the movement of the UAV during the flight also caused the angular coverage per pixel to be slightly different for each spectral band ( Figure 5). For the fitting of a BRDF model, this does not have to be a problem as long as a large amount of observation angles are available per pixel.
To capture a representative amount of potato rows and spaces between the potato rows in each ground pixel, we resampled the original ground pixels (~8 cm) to 5 m ground pixels. For more homogeneous surfaces or surfaces with smaller height differences, a smaller or even the original ground pixel size can be used. The proposed methodology in this study can easily be applied to-which is more commonly available-multispectral camera systems. We used the Rikola camera system because that was the camera we had available. Multispectral cameras, which in general have a higher spatial resolution, allow for analysis at a higher detail level. If local legislation allows it, these cameras can be flown at greater heights above ground level to cover a larger surface area, while maintaining a high spatial resolution. Multispectral camera systems with bands in both the visible and NIR capture the most pronounced changes in anisotropy that relate to increased canopy cover in a potato crop, as can be seen in the Θ and ρ 0 parameter spectra in Figure 9 and the correlation spectra in Figure 11.
We studied the magnitude of anisotropic reflectance effects by analyzing anisotropy parameters obtained by fitting the RPV model through the multi-angular observations of each georeferenced ground pixel. Observed trends were a backward scattering anisotropy for the whole study area, which was most pronounced in the visible wavelength region and less pronounced in the NIR region, which is typical for vegetation targets [2,3].
In addition, we studied the effects of canopy cover and LAI on reflectance anisotropy. On the first measuring day, we observed an increase in backward scattering intensity with increased canopy cover at visible wavelengths. This might be related to the higher amount of gaps in the potato rows on this day, which in turn resulted in weaker shadows and thus weaker backward scattering intensity. A higher canopy cover on this day occurred at plots where the potato rows had fewer gaps, which resulted in stronger shadow effects. On the second measuring day, when there were fewer gaps between rows, we observed a decrease in backward scattering intensity with an increase in canopy cover. This decrease is likely the result of the filling up of space between consecutive rows due to potato plants growing into each other (the potato plants had fallen over due to long and unstable haulms in combination with unfavorable weather conditions), which resulted in a decrease in shadows created by the rows. A decrease in backward scattering intensity with canopy closure was observed in a previous UAV-based anisotropy study of potato canopy development [36].
This study demonstrates that pixel-wise multi-angular observations can be obtained from frame-based cameras by calculation of the VAAs and VZAs of georeferenced ground pixels that were captured from multiple camera positions. The information contained in these multi-angular observations is valuable for correction of anisotropic reflectance effects, but also provides additional information to spectral data on vegetation characteristics. Moreover, knowing the observation geometry of every pixel makes it possible to select and analyze only observations that were taken from, for example, the backward scattering direction. In this direction, observations have been demonstrated to be more useful for LAI prediction [56]. Besides this, in the backward scattering direction, shadows are minimized and observations are more sensitive to vegetation pigments [57], which makes frame-based measurements in combination with their observation geometry also interesting for the study of biochemical vegetation parameters.

Conclusions
The overlap of images that are collected during a typical UAV mapping flight using a frame camera provides multi-angular views of georeferenced ground pixels. In this paper, we describe how to extract the observation geometry of these multi-angular views. The results of this paper Remote Sens. 2017, 9, 417 20 of 23 demonstrate that, in clear sky illumination conditions, strong anisotropy effects occur in measurements with off-nadir viewing when using a frame camera.
By fitting the RPV model through the multi-angular measurements, we parameterized and interpreted the anisotropic reflectance effects captured during two UAV flights in the growing season of 2016 for a potato field that contained eight experimental plots with different LAI and canopy cover levels. Analysis of the anisotropy patterns indicated that in general the potato crop had a higher reflectance in the backward scattering direction. An increase in canopy cover during the growing season resulted in a reduction of this backward scattering, which can be attributed to the diminishing row structure with increasing canopy cover, and thereby weaker row-induced shadow effects.
The results of this study indicate that the parameters describing the anisotropy with the RPV model contain information on structural characteristics of the potato crop such as LAI and canopy cover. This suggests that including such anisotropy information may lead to more accurate estimates of these characteristics. Future research will focus on how to use this anisotropy information to better estimate vegetation parameters, for example, by using radiative transfer models. Moreover, since the method described in this paper can in theory be used to study the anisotropy of any surface, we will apply it to other vegetation, crop, and soil targets as well. Special focus will be on the influence of canopy development on the anisotropy signal.