AUTOMATIC EXTRACTION OF SOLAR AND SENSOR IMAGING GEOMETRY FROM UAV-BORNE PUSH-BROOM HYPERSPECTRAL CAMERA

: Calculating solar-sensor zenith and azimuth angles for hyperspectral images collected by UAVs are important in terms of conducting bi-directional reflectance function (BRDF) correction or radiative transfer modeling-based applications in remote sensing. These applications are even more necessary to perform high-throughput phenotyping and precision agriculture tasks. This study demonstrates an automated Python framework that can calculate the solar-sensor zenith and azimuth angles for a push-broom hyperspectral camera equipped in a UAV. First, the hyperspectral images were radiometrically and geometrically corrected. Second, the high-precision Global Navigation Satellite System (GNSS) and Inertial Measurement Unit (IMU) data for the flight path was extracted and corresponding UAV points for each pixel were identified. Finally, the angles were calculated using spherical trigonometry and linear algebra. The results show that the solar zenith angle (SZA) and solar azimuth angle (SAA) calculated by our method provided higher precision angular values compared to other available tools. The viewing zenith angle (VZA) was lower near the flight path and higher near the edge of the images. The viewing azimuth angle (VAA) pattern showed higher values to the left and lower values to the right side of the flight line. The methods described in this study is easily reproducible to other study areas and applications.


INTRODUCTION
Remote sensing has proved to be highly effective and efficient in studying a diverse variety of natural and ecological resources. Other than satellite and aerial remote sensing, recent advances in Unmanned Aerial Vehicles (UAV) and sensor technology has opened more opportunities to study vegetation dynamics, specifically in agricultural applications (Maddikunta et al., 2021). Since UAVs can be flown at lower altitudes than satellites or aircrafts, the resulting products offer higher spatial resolution and with more accurate canopy spectra (Tao et al., 2020). The canopy spectra can be used to model or represent different plant traits. For instance, different vegetation indices (e.g., normalized difference vegetation index, NDVI) can indicate overall crop health that improves precision agriculture practices (Radoglou-Grammatikis et al., 2020). Additionally, UAV sensors offer highthroughput plant phenotyping that accelerates current crop breeding operations (Song et al., 2021). Moreover, UAV-based imageries can be used to train advanced machine learning models, which predict various crop traits, disease, yield, and seed quality at plot-level (Bhadra et al., 2020;Maimaitijiang et al., 2020;Nguyen et al., 2021).
Hyperspectral sensors can collect reflected spectra from crop canopies with higher spectral resolution. A typical hyperspectral image (HSI) often contains hundreds or even thousands of bands for a wide range of wavelengths. Generally, the wavelengths can vary from Very Near Infrared (VNIR, 400-1000 nm) to Shortwave Infrared (SWIR, 900-2500 nm) with different spectral resolution (1 to 10 nm). Plants reflect electromagnetic radiation, which contains information about their biophysical composition and physiological status (Segarra et al., 2020). Numerous studies * Corresponding author have utilized the broader range of HSI products to study different characteristics of plants and vegetation (Mariotto et al., 2013;Fernandes et al., 2015;Banerjee et al., 2020;Wang et al., 2021).
The quality of HSI-based inference heavily depends on the accuracy of HSI post-processing techniques. Generally, the HSI sensor provides the raw Digital Number (DN) or radiance (in W•sr -2 •m -2 ), which is then converted to unitless top-ofatmosphere (TOA) reflectance and surface reflectance (SR). Empirical Line Method (ELM) is the widely used calibration technique to convert radiance into SR using different calibration targets on the ground (Markelin et al., 2008;Wang and Myint 2015;Ortiz et al., 2017). The principal assumption behind this technique is that the objects on the ground represent a Lambertian surface, which appears uniformly bright from all directions of view and reflects the entire incident light (Mao et al., 2020). However, the crop canopy architecture is far from being a Lambertian surface and exhibits anisotropic effects (Jiao et al., 2014). Therefore, several studies have identified that multiple viewing angles or viewing geometry of sensors play an important role in the pixel-level SR (Vermote et al., 2009;Zhang et al., 2014). For example, Galvao et al., (2009) retrieved highly accurate Vegetation Indices (VIs) from Hyperion and MODIS satellite data when using backward observations. Similarly, Gu et al., (2015) found improved Leaf Area Index (LAI) estimation accuracy from backward observations compared to forward observations in CHRIS/PROBA data. Huang et al., (2011) demonstrated that multi-angular hyperspectral observations could retrieve the vertical distribution of chlorophyll content in winter wheat.
In terms of UAV-based observations, several studies have utilized snapshot (or frame) hyperspectral cameras to derive multi-angular spectral information. For instance, Roosjen et al., (2018) achieved improved results in estimating LAI and chlorophyll content of potato by using multi-angular data. They introduced a goniometer-based simulation method for HSI footprints with high overlaps. The multiple viewing angles were converted to zenith and azimuth angles, which were used to simulate PROSAIL spectra and derive better LAI and chlorophyll retrieval accuracy. Similarly, Mao et al., (2020) found that the effect of multi-angular observations was significant in deriving VIs for soybean and maize. They also extrapolated different viewing angles from a snapshot hyperspectral camera mounted with a UAV and corrected for the Bi-directional Reflectance Function (BRDF) effect. Therefore, the availability of multiple solar-sensor zenith and azimuth angles is highly important to accurately study different plan characteristics.
Alternative to snapshot cameras, push-broom hyperspectral sensors (or line-scanner sensors) are now widely used with UAVs. The push-broom sensor captures one line per exposure that forms one image line after the other (Barreto et al., 2019). Therefore, push-broom sensors can outperform snapshot cameras, as the latter systems require a compromise between spatial coverage, spatial resolution, and spectral resolution (Aasen et al., 2015;Yi et al., 2021). However, extracting the solar-sensor zenith and azimuth angles from a push-broom sensor is not as straightforward as snapshot cameras. While the snapshot camera provides 2D scenes captured across overlapping flight lines in relatively higher time interval, the push-broom sensor captures line by line 1D spectra across its flight path. Due to the line-by-line scanning mechanism, push-broom sensors suffer from wind-related motions during data acquisition (Jaud et al., 2018). As a result, push-broom sensors require high accuracy Global Navigation Satellite System (GNSS) and Inertial Measurement Unit (IMU) onboard the UAV to ortho-rectify the lines and generate a geometrically accurate hyperspectral cube (Yuan and Zhang 2008). Due to the availability of GNSS/IMU system onboard the platform, the solar-sensor geometry can be directly calculated using linear algebra and spherical trigonometry. Therefore, the objective of this study is to develop an automated framework that can calculate the solar-sensor zenith and azimuth angles for each pixel in a hyperspectral cube collected by a push-broom UAV scanner with cross-grid flight pattern.

Experimental Setup
The experiment was setup in the Planthaven Farms at OFallon, Missouri, United States ( Figure 1). The site was located slightly northeast from Saint Louis city close to the Mississippi River to the north. The field was planted with 220 rows of maize on May 25, 2021, where 2 rows were marked as one plot. Total 55 different genotypes or cultivars of maize were planted with 2 replicas. The field was approximately 75 m long and 20 m wide. During the growing season, average temperature was between 23-24°C and average annual precipitation was 1092.2 mm for the study area.

UAV Flight
A DJI M600 Pro UAV was used to collect the hyperspectral data for the study area ( Figure 1c). The UAV was equipped with a Headwall Nano-Hyperspec VNIR push broom camera (Headwall Photonics, Massachusetts, United States), a FLIR Vue Pro thermal camera (FLIR Systems, Oregon, United States), and an APX-15 GNSS/IMU (Applanix Corporation, Ontario, Canada) unit all attached to a DJI Gimbal (Figure 1d). The APX-15 UAV GNSS/IMU records the precise time, position, and orientation of the sensor at 200 Hz interval. Full specifications of the sensor and the GNSS/IMU unit is provided in Table 1. Two UAV flights were conducted on July 20th and August 4th of 2021. Each flight was planned in a cross-grid pattern (Figure 1f) in UgCS mission planning software (v4.0.187, SPH Engineering, Latvia) with 4 length wise and 9 width wise lines, resulting in total 13 hyperspectral cubes. The altitude and velocity for both flights were 50 m and 3 m/s. The ground sampling distance (GSD) was found 3.01 cm from both flights.

METHODS
The extraction of solar-sensor geomtery contain three major parts ( Figure 2): 1) hyperspectral cube processing, 2) locating viewing point for each pixel, and 3) calculating solar-sensor geometry. The process was automated using Python libraries which are available in a public repository with test datasets (https://github.com/remotesensinglab/uav-solar-sensor-angle).

Figure 2.
Overview of methods for the extraction of solar-sensor zenith and azimuth angles, i.e., (a) hyperspectral cube processing involves converting digital number (DN) to radiance, reflectance, and ortho-rectified images, (b) locating corresponding sensor location for each pixel in a hyperspectral cube, and (c) calculating three solar-sensor zenith and azimuth angles.

Hyperspectral Cube Processing
The raw hyperspectral data cubes (HSI cubes) contain all the data as digital number (DN) in the raw images. The DN values does not provide any meaningful information about the scene, rather provides 12-bit numeric values collected by the sensor. Each flight lines resulted in a separate HSI cube. Before the flight, a dark reference data cube was collected and stored within the sensor. The DN values were converted to radiance (W•sr -1 •m -2 ) values using the dark reference and the vendor provided software named Headwall SpectralView software (v3.1.4, Headwall Photonics, Massachusetts, United States). The next step was to convert the radiance to reflectance values (%) by using a 56% reflectance tarp used in the field during the data collection. Finally, the cubes were ortho-rectified using the same software suit and APX-15 high-resolution GNSS/IMU dataset. The GNSS/IMU dataset was corrected using a nearby Continuously Operating Reference Stations (CORS) base station located in OFallon, Missouri and the vendor provided POSPAC UAV software suit (v8.7, Applanix Corporation, Ontario, Canada). The overall procedure is illustrated in Figure 2a.

Locating Viewing Point for Each Pixel
The viewing point in terms of each pixel was required to calculate both sensor zenith and azimuth angles (Figure 2b). First, the GNSS data was extracted from APX-15 device and converted to an ASCII text file which contained latitude, longitude, and timestamp information. Additionally, the coordinates of each pixel was calculated by converting the raster data into a geospatial text file using GDAL v3.3.1 (GDAL/OGR 2020). The raster image also had longitude, latitude, and timestamp information. Therefore, only the corresponding GNSS observations for a HSI cube were filtered by matching the timestamp from the cube and GNSS points. Finally, the GNSS points which had the shortest distance from each pixel location were identified by representing the point pairs in a matrix form. It was done by calculating Euclidean distance from each pair of pixels and GNSS coordinates. The information was preserved in a text file as comma-separated value (CSV) format, which contained the unique ID of the closest GNSS point for every pixel in the HSI cube.

Solar-Sensor Angle Calculation
Overview of solar zenith angle or viewing zenith angle (VZA, θ V ), solar zenith angle (SZA, ), sensor or viewing azimuth angle (VAA, ϕ V ), and solar azimuth angle (SAA, ϕ S ) calculations in terms of cross-grid UAV with push-broom hyperspectral sensor are illustrated in Figure 2c.

Solar Zenith Angle (SZA):
The solar zenith angle (SZA) is similar to the VZA, but instead of the sensor as the moving vector, the position of the sun becomes the point of interest. SZA is a function of the raster location coordinates (longitude and latitude) and time of the day. SZA (θ S ) can be calculated from Solar Elevation Angle (α S ) using Equation 3.

= 90 −
(3) sin = sin sin + cos cos cos ℎ where, ϕ is the latitude of the location, δ is the solar declination angle and h is hour angle. Declination angle (δ) is the angle between the line joining the centers of the Sun and the Earth and its projection on the equatorial plane.
The pixel coordinates were attached with corresponding sensor points and each sensor point included the time information in UTC format. Also, the coordinates were converted from UTM to a geographic coordinate system (World Geographic System 1984), so the values were available as latitude and longitude. A python package called PVLIB (v0.9.0) was used to calculate δ, h and eventually θ S for all pixel coordinates.

Solar
where ϕ, δ, h and α S are the latitude, solar declination angle, hour angle and solar elevation angle, respectively.

Viewing Zenith Angle (VZA):
The viewing zenith angle (VZA) is the angle between the vector from sensor and raster point (VR ⃗⃗⃗⃗⃗⃗ ), and the surface normal (Z ⃗ ⃗ ) from the raster point (also known as zenith), which can be defined as θ V . The UAV was flown at a 50 m altitude for the whole mission. Therefore, a perpendicular vector from the sensor point to the XY surface can be drawn as VV ⃗⃗⃗⃗⃗⃗ , where VV ⃗⃗⃗⃗⃗⃗ is 50 m. The angle between VR ⃗⃗⃗⃗⃗⃗ and RV ⃗⃗⃗⃗⃗⃗ is known as Viewing Elevation Angle (α V ) and can be calculated using Equation 1. (1) where the coordinates of R and V are (x r ,y r ) and (x v ,y v ), respectively, calculated in a Universal Transverse Mercator (UTM) projection system. Therefore, RV ⃗⃗⃗⃗⃗⃗ can be calculated as the Euclidean distance between R and V . If α V is known, then θ V can be calculated using Equation 2.

= 90 − (2)
For every pixel coordinate, corresponding θ V values were calculated and converted to degrees.

Viewing Azimuth Angle (VAA):
Azimuth angle can be calculated in the XY plane, where it is the clockwise angle between a point of interest and the true north (Y ⃗⃗ ). For calculating the sensor or viewing azimuth angle (VAA, ϕ V ), an arbitrary north vector for any pixel coordinate (x r ,y r ) was created by adding 100 m to y r (Figure 2c2). The VAA can be calculated using Equation 7.
where a is the true north vector and b ⃗ is the vector between the raster point and corresponding sensor point.

RESULTS AND DISCUSSIONS
The SZA (θ S ), SAA (ϕ S ), VZA (θ V ), and VAA (ϕ V ) were calculated for all 13 HSI cubes, but we will discuss only the 1 st HSI cube (Figure 3). Additionally, the descriptive statistics of the angles are provided in Table 2.   Table 2. Descriptive statistics of the angles for the 1 st HSI cube.

Solar Angles
The SZA and SAA shows different angular pattern in the resulting rasters (Figure 3b and 3c). The SZA started decreasing along with the flight direction, whereas the SAA started to increase along with the flight direction. The total duration for capturing this HSI cube was 42.32 seconds, which resulted in low standard deviation for the angular values of SZA and SAA ( Table  2).
The verification of our calculation of solar angles was done by calculating the solar position based on National Oceanic and Atmospheric Administration (NOAA) Solar Calculator (NOAA 2021), which is an online tool to calculate approximate solar position in terms of coordinates and local time. To verify our result with the NOAA Solar Calculator, 5 randomly selected points were selected, and corresponding solar angles were extracted. Table 3 shows the SZA and SAA calculated based on NOAA Solar Calculator and our method, and the absolute differences observed.  Table 3. Comparison of solar angles between NOAA Solar Calculator and our method. Abs. Diff. indicates the absolute differences between two methods and T is the order of points mentioned in Figure 3.
The SZA and SAA values calculated by NOAA Solar Calculator and our method showed slight differences at the decimal level. Since we used highly accurate PVLIB Python library to calculate the solar angles, we could provide coordinates and time information up to any decimal level possible. For instance, the time information in our method had 6 decimal places for second values. On the other hand, the NOAA calculator could only take the second values as integer. Moreover, NOAA (2021) indicates that due to the variations in the atmospheric conditions and uncertainty in the algorithms, there could be slight differences in the solar position calculations. These could be attributed to the slight differences in the solar angle values. However, having precise coordinate and time information in the angle calculation is highly preferable for remote sensing applications, specifically in HSI-based processing.

Sensor Angles
The pattern of VZA (Figure 3d) can be explainable in terms of the flight path. Since the flight path runs through the middle of the cube, the VZA values are close to zero near the flight path and starts increasing at the edge of the image. Since this is zenith angle from the sensor, there should be higher angles at the edge rather than the middle.
The VAA shows comparatively larger range of angular values (Figure 3e). Since azimuth angle is calculated as the clockwise angle from the north vector of each raster pixel, the right side of the cube resulted with smaller angular values, whereas the left side comprised of higher values. However, the pattern in the VAA raster may seem binary, but the inset map on Figure 3e shows an enlarged portion of the right side. The inset map shows that there exists angular variation along the flight path and the variation can be seen perpendicular to the flight line. Therefore, the standard deviation is the highest for VAA with larger range of angular values (Table 2).

Limitations
The major issue encountered in this study was the lack of camera calibration. We used the vendor provided software (Headwall Spectral View v3.1.4) to ortho-rectify the HSI cubes. However, when the cubes were plotted in a GIS environment, it was noted that the overlapping regions from two consecutive HSI cubes did not exactly match. Therefore, the HSI cubes were georeferenced with a Light Detection and Ranging (LiDAR)-derived RGB point cloud using 6 control points. The LiDAR mission was also flown on the same days the HSI missions were performed. The LiDAR UAV point cloud was corrected using a GNSS base station established during the data collection time and the vendor provided software named, Phoenix LiDARMill (v2.0, Phonix LiDAR Systems, Texas, United States). After correction, the position accuracy was around ±0.1 cm. When the HSI cubes were georeferenced with the RGB point cloud, corresponding cubes matched properly with each other.
However, the problem can be solved by performing a camera calibration. Probably the internal operating parameters (IOPs), boresight angles or the lever-arm offsets were changed from the initially approximated values by the vendor. LaForest et al., (2019) performed similar camera calibration technique to perform time-delay adjustment to similar type of HSI UAV platform that had accurate GNSS/IMU information. The same methodology can be applied to our study to improve upon the ortho-rectification of the HSI cubes as well. Therefore, careful considerations should be made when working with push-broom sensors and calibration flights should be conducted using randomly placed ground control points (GCPs) on the ground.

CONCLUSIONS
The study demonstrates a simple methodology for calculating solar-sensor zenith and azimuth angles for a push-broom hyperspectral sensor equipped in a UAV. The results show that the method can deliver all the angles in raster format, which can be very helpful to perform BRDF corrections or radiative transfer model-based applications in remote sensing of vegetation. If this work is needed to be reproduced for other study areas using similar sensor and platform, then it will be easy to do so by utilizing the automatic Python workflow developed from this work. In future, we will improve the camera calibration issue incurred in this study and apply these angle rasters to perform radiative transfer modeling-based applications for plant phenotyping.