Data-Driven Artificial Intelligence for Calibration of Hyperspectral Big Data

Near-earth hyperspectral big data present both huge opportunities and challenges for spurring developments in agriculture and high-throughput plant phenotyping and breeding. In this article, we present data-driven approaches to address the calibration challenges for utilizing near-earth hyperspectral data for agriculture. A data-driven, fully automated calibration workflow that includes a suite of robust algorithms for radiometric calibration, bidirectional reflectance distribution function (BRDF) correction and reflectance normalization, soil and shadow masking, and image quality assessments was developed. An empirical method that utilizes predetermined models between camera photon counts (digital numbers) and downwelling irradiance measurements for each spectral band was established to perform radiometric calibration. A kernel-driven semiempirical BRDF correction method based on the Ross Thick-Li Sparse (RTLS) model was used to normalize the data for both changes in solar elevation and sensor view angle differences attributed to pixel location within the field of view. Following rigorous radiometric and BRDF corrections, novel rule-based methods were developed to conduct automatic soil removal; and a newly proposed approach was used for image quality assessment; additionally, shadow masking and plot-level feature extraction were carried out. Our results show that the automated calibration, processing, storage, and analysis pipeline developed in this work can effectively handle massive amounts of hyperspectral data and address the urgent challenges related to the production of sustainable bioenergy and food crops, targeting methods to accelerate plant breeding for improving yield and biomass traits.


I. INTRODUCTION
R EMOTE sensing imagery in the visible, near-infrared (VNIR) and shortwave infrared (SWIR) (400-2500 μm) domains are affected by noise and uncertainty attributed to the atmosphere (e.g., effects due to molecular scattering and absorption), solar and sensor viewing geometries, terrain, environmental context, and mechanical degradation of sensor optics over time. There is no unified approach to address these challenges as they are variable over timescales and geographies. For example, atmospheric constituents such as haze, aerosols, clouds, and water vapor are all different over time and space. Therefore, atmospheric correction routines should be versatile and adjustable to account for these differences. Additionally, massive amounts of structured and unstructured datasets from various sensors (multispectral, hyperspectral, thermal, SAR) and platforms (e.g., drones, polar or solar orbiting, and geostationary satellites) have accumulated over decades, and the remote sensing community has been drowning in big data. These diverse datasets are crucial for various applications and informed decision-making. However, they may carry nonuniform spectral information that can be different both in magnitude and shape for the same target partly due to the differences in spectral sensitivity, signalto-noise ratio (SNR), and ground sampling distance (GSD) of these sensors but most importantly disturbances resulted from atmospheric scattering and absorption and anisotropy associated with terrain heterogeneity. If not accounted for, bias from local atmospheric and terrain conditions severely limits the use of these datasets for spatiotemporal analysis, quantitative estimation of biophysical and chemical properties of matter, and detection of change in global studies [1]. Methods developed for one site or data source may not be This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ applicable to a different site or datasets. This emphasizes the need for the development of robust, artificial intelligence (AI)-powered atmospheric and radiometric processing (i.e., conversion of sensor radiance to surface reflectance) tools that allow datasets from different sensors and time to demonstrate the same spectral features (e.g., absorption or reflectance features in surface reflectance) in which a target is known. Beyond atmospheric corrections, other parameters such as inference from background soil and shadows need to be addressed to derive accurate and reliable spectral measurements.
Since the inception of imaging spectroscopies, issues related to accurate and reproducible data have been a primary concern due to the impact of illumination and atmospheric effects on the resulting spectral measurements. The process is twofold: first, top-of-atmosphere (at sensor) radiance is calculated, and atmospheric properties (e.g., transmittance, absorption, and scattering parameters) are estimated along with surface reflectance spectra reflective of surface physical, chemical, or compositional properties. To reduce these effects from spectroscopy data, several methods have been proposed. Early radiometric calibration methods such as the empirical line method (ELM) were developed for the correction of aerial and satellite-based multispectral and hyperspectral data to transform raw digital numbers (DNs) or radiance to at-surface reflectance [2], [3]. ELM provides a simple solution to correct imagery data by establishing an empirical relationship between two or more ground targets with a known reflectance factor and the observed radiance values. This reduces the calibration for each spectral band to an uncomplicated regression equation. The result is the calculation of gain and offset coefficients that can be applied to all surfaces within the hyperspectral data, assuming uniform atmospheric conditions. The gain characterizes the sensor response of reflectance factor per unit radiance and the offset characterizes the sky path radiance between the sensor and the surface for a sensor system calibrated to radiance [4]. Although the ELM is by far the most common calibration procedure, several problems have been identified [5]. One of the major flaws is the assumption of a Lambertian surface, which is rarely the case in nature. Another flaw is that at-sensor radiance is not easily obtained for airborne or satellite imagery since the relationship between surface reflectance and DNs is not always linear and may change across bands.
To estimate a more complete radiometric calibration, methods based on radiative transfer simulations have been developed. These methods calculate atmospheric effects based on the physical radiative-transfer model and the atmospheric conditions when the hyperspectral image is acquired by the sensor. Measured surface and atmospheric properties from the target location can then be used in a radiative transfer code such as MODTRAN [6]- [8] and 6S [9] to predict the band-averaged, top-of-atmosphere radiance at each spectral band. Although very accurate, radiometric simulations are time-consuming and not scalable computationally in an automated workflow.
Herein lies the opportunity for AI and data-driven approaches to help understand and model the world, by extracting patterns from big data. Data-driven approaches have been deployed in a variety of fields where data is ubiquitous, and manual annotation or modeling is difficult. For example, data-driven approaches have yielded recent successes in drug discovery [10]- [12], health care [13]- [15], and urban planning [16], [17]. The AI approaches and algorithms deployed in different problem domains vary, but the general idea is to use both supervised (labeled) and unsupervised (unlabeled) approaches to gaining insights from large datasets to understand and model complex systems. Such approaches have been deployed for hyperspectral data analysis pipelines. For example, recently, neural network approaches [18], [19] and optimal estimation (OE) [1], [20], [21] have gained much attention in the atmospheric correction of airborne and satellite imagery, which enables the development of a generalized approach for various surface types or biome (aquatic, terrestrial, or/and inland water surfaces) [21].
In this article, we develop an AI-based pipeline for organizing and calibrating super high-resolution hyperspectral big data captured for the purpose of improved bioenergy crop phenotyping. In 2015, the U.S. Department of Energy funded a program for Transportation Energy Resources from Renewable Agriculture (TERRA) through the Advanced Research Projects Agency-Energy (ARPA-e) to address the urgent challenges related to the production of sustainable bioenergy crops, targeting methods to accelerate plant breeding for improved yield and biomass traits through advanced technologies for high-throughput crop phenotyping. Given that food and bioenergy crops are improved to be more waterand nutrient-use efficient through plant breeding methods which require accurate measurements of phenotypic traits and associated genotypic traits, understanding and predictability of complex traits like yield and biomass can be improved with increased throughput, accuracy, and availability of field phenotyping data [22], [23]. As part of the project, the TERRA-REF Field Phenotyping Reference Platform (aka the field scanner) was constructed at the University of Arizona Maricopa Agricultural Center and operated during a total of nine field trials between 2016 and 2019 to generate publicly available plant phenotype data on diverse sets of sorghum and durum wheat germplasm. The gantry "camera box" supports five primary imaging systems: a 3-D laser, thermal infrared, stereo RGB cameras, a chlorophyll fluorescence camera, and a hyperspectral imaging system. The hyperspectral datasets generated by two Headwall Photonics instrument systems are valuable reference data since high-resolution hyperspectral imaging systems are a promising technology for monitoring and evaluating the physiological status of plants and plant systems across scales from microscopic to ecosystem levels [24]. The field scanner is an automated field "robot" that can be preprogrammed to collect data and is capable of producing terabytes of data per day, making calibration of hyperspectral big data an enormous and important challenge. While the importance of collecting plant trait data in field environments is recognized given that plant growth and yield are the result of complex interactions with the dynamic environment, the same features of a dynamic and diverse range of environmental variables that increase the importance of field-collected data present immense challenges to data collection, calibration, processing, and interpretation. Manual handling of these massive collections of data is impossible. It is crucial to develop data-driven and AI-powered workflows and algorithms that are capable of quality checking and early warning of data quality, versatile and generalized calibration pipeline for radiometric and atmospheric correction, soil/shadow removal, and bidirectional reflectance distribution function (BRDF) normalization.
The objectives of the work presented in this article are to: 1) develop a suite of algorithms for calibration of gantry hyperspectral data including radiometric calibration, BRDF correction and reflectance normalization, and soil/shadow masking and 2) describe a data-driven AI for automation of a calibration pipeline for hyperspectral big data.
It is worth noting that this article represents a novel concept of AI for remote sensing field. AI is embedded in every step in the processing chain including data collection, transfer and management, calibration, and phenotyping with various algorithms as opposed to a single AI algorithm. For example, a fully automated field robot collects data in a fully automated imaging procedure in which aperture opening, platform speed, and the location of scan are determined by real-time atmospheric conditions and plant growth. The data are then transferred to the cloud automatically where image quality assessment, calibration, and feature extraction take place. Image quality assessment algorithms inspect image quality and usability and trigger a warning message to the database manager if there are any quality issues that are fully based on AI. Then the calibration part utilizes machine learning (which is a subfield of AI) for almost all processing steps.

A. Study Site and the Gantry Hyperspectral System
The study area ( Fig. 1) is the TERRA-REF field scanner [ Fig. 1 (b)] field site at the Maricopa Agricultural Center, Maricopa, AZ, USA (33.070 • N, 111.974 • W, elevation 360 m). The scanner is constructed on parallel rails spaced 28 m apart (providing a scannable field width of 23.8 m) over a length of 218 m in the original design, expanded another 170 m in length in 2018, to enable 0.47 hectares (0.84 hectares after rail expansion) of field plots to be scanned noninvasively from approximately 2.5 m above the plant canopy. Field plots in the scannable area can be planted, harvested, and maintained using full-size field equipment. The original rail length allowed a total of 700 field plots to be scanned (two-row plots, each 3.5 m in length). Camera box motion during scans is programmed in three dimensions (along the short and long axes of the field as well as height above ground), and scans are semiautomated requiring minimal operator intervention. The velocity of the camera box during a scan and the proportion of plot coverage are programmable variables such that scans across the field vary in the amount of time to cover the field. Most scans take hours to complete and span a wide range of diurnal solar angles and cloud conditions. Scans can be set to run throughout the crop cycle from planting to harvest and during most weather conditions with the exception of imminent risk of lightning or high wind speeds.

B. Gantry Hyperspectral System
The hyperspectral imaging system in the field scanner is comprised of two HyperSpec Inspector line-scanning cameras (Headwall Photonics) in environmental enclosures mounted in the camera box and two associated downwelling spectroradiometers (Ocean Optics) mounted on an arm from the top of the gantry support. The VNIR sensor is a scanning imager for the visible and near-infrared range from 389 to 1000 nm, while the SWIR scanning imager covers the range from 895 to 2501 nm. During the data collection time span over multiple field seasons, the VNIR and SWIR scanners were returned to Headwall Photonics for upkeep and calibration, resulting in slight modifications to the sensors and data outputs. Both scanning imagers have a slit size of 25 μm. The VNIR camera has a spectral resolution of 0.64 nm, 955 (before August 2018)/939 (after August 2019) wavebands, spectral rows 2160 pixels, nominal pixel pitch 0.0065 mm, and average dispersion of 98.4 nm/mm. The SWIR camera has a spectral resolution of 5.9 nm, 273 (before August 2018)/275 (after August 2018) wavebands, spectral rows 288 pixels, nominal pixel pitch 0.024 mm, and average dispersion of 246.4 nm/mm. The reflected light measured in photons (aka DN) by the hyperspectral scanning imagers needs to be further processed, which is a process called radiometric calibration that converts photons to radiances in watts per square meter per steradian (W·sr −1 ·m −2 ·um −1 ) or unit-less reflectance factor. The common approach for radiometric calibration is to use a known target reflectance by a ratio of upwelling radiances (DN) to downwelling irradiances measured by a spectroradiometer pointed upward. During the time from April 2016 to February 2019, a single spectroradiometer was in operation to measure the spectral range 350-800 nm at an optical resolution of 3 nm (Ocean Optics, STS-VIS-L-50-400-SMA), recording one measurement every 5 s (0.2 Hz). In early 2019, when funding became available, the single spectroradiometer was replaced with two instruments to cover the combined spectral ranges of 350-1000 nm, 3648 pixels (Ocean Optics FLAME-T-VIS-NIR-ES), and 900-2500 at a spectral resolution of ∼6.3 nm (full width at half maximum, FWHM) (Ocean Optics NIRQUEST512-2.5). The field scanner datasets generated by the hyperspectral imaging system during frequent scans and the spectroradiometers during continuous logging provided the raw datasets that subsequently need to be georeferenced, linked with associated metadata, quality checked, and calibrated before they can be fully available as an open data resource for users, including the public and phenotyping community.
The TERRA-REF field scanner and imaging instruments were constructed to be able to collect data in diverse and often harsh environmental conditions, providing great opportunities to monitor plants in stressful conditions, but also requiring continual system operational monitoring and adjustment with occasional repairs. Scans have been conducted during summer day temperatures of 48 • C in June 2017 and during winter night temperatures of −3.9 • C in January 2019. The inevitable adjustments and repairs that are necessary to optimize instrument function across a wide range of conditions have a downstream impact on data processing. For example, processing algorithms that work well for a specific camera exposure setting may not work across all exposure settings needed to cover conditions from clouds to full sun and across diurnal sun angles. Data outputs from the hyperspectral imaging system are sensitive to light intensity, wind, and variability in plant canopy structure, all of which intersect to create complex mixes of sunlit and shaded plant surfaces. In summary, while data from controlled-environment studies allow more uniform environmental conditions for easier data collection, processing, and interpretation, the fact that plant physiological efficiency is an integration of abilities to respond to dynamic environments establishes the need for datasets from field trials [25].

A. Calibration Workflow
Hyperspectral data processing mainly includes radiometric calibration, BRDF correction, shadow masking, soil pixel removal, georeferencing, NetCDF format conversion, and plotlevel feature extraction modules. Different processing levels were categorized as level 0, level 1, level 2, and level 3 (Fig. 2). Headwall hyperspectral sensors integrated on the field scanner collected high-resolution VNIR and SWIR imagery in digital counts/numbers [DN, which is similar to a photon counter modulated by the spectral response function (SRF) of the specific sensor] and were saved as an HDR file, which is the raw imagery data in the hyperspectral data processing workflow.
First, radiometric calibration models were established to convert the VNIR and SWIR imagery raw values to a hyperspectral reflectance factor using synchronously recorded irradiance data from real-time downwelling irradiance sensors mounted on the top of the field scanner. Second, BRDF correction models were built to normalize the reflectance deviation due to various illumination conditions and sensor viewing angles during data collection. A machine learning-based Overall workflow of hyperspectral imagery processing pipeline. Level 0 raw imagery is automatically uploaded to a cloud storage in real-time, then DNs are converted to a reflectance factor (Level 1) by empirical correction methods. Level 1 reflectance factor data are then normalized to minimize the effects of illumination (changes in sun elevation) and viewing geometry (variations of measurements within the field of view along the scan width) producing Level 2 analysis-ready reflectance factor data without the soil and shadow component. Level 3 data are plot level features including various spectral indices and plant phenotypic traits. shadow-masking model was developed to exclude the shadow cast by the field scanner's frame and camera box. A soil mask was generated to exclude soil pixels from further processing. Once completed, the geospatial information was assigned, and the high-quality reflectance imagery data was converted to a Network Common Data Form (NetCDF) format. Finally, plot-level reflectance imagery was extracted, along with a variety of spectral features such as vegetation indices that were computed at the plot level. The details of each processing module will be described and explained in the following sections.

B. Radiometric Calibration
Radiometric calibration transforms the source data with digital counts/DNs or the number of photons measured by hyperspectral spectroscopy to radiance or reflectance factor values that has a physical meaning, which is often used to conduct a quantitative analysis of the spectral, spatial, and temporal characteristics of the surface targets [26]. The radiometrically calibrated hyperspectral reflectance factor data allows an efficient way to derive plant canopy spectral information/features such as vegetation indices, which are important indicators of plant growth status, and representations of plant biochemical and biophysical characteristics that contribute to yield and can indicate plant stress [27]- [31].
Dividing the hyperspectral sensor captured radiance value by irradiance value is often used for radiometric calibration to produce reflectance factor data. Radiance values could be converted from the hyperspectral sensor source data in DNs or intensity values, which requires sensor calibration files, including SRFs, sensor gain, and offset information derived in a laboratory setting. However, such information is often sensor specific, and proprietary to the vendor is not readily available in our case. Reflectance factor could also be obtained using an ELM [26] established between reference targets with known reflectance factor values and photon counts (also known as DNs) values measured by an imaging sensor. In our case, the hyperspectral sensor calibration parameters that allow converting DNs to radiance were not available. It should be mentioned that the diurnal or temporal changes of illumination conditions during data collection make it impossible to use in situ reference targets and automate the process because it is not practical to reliably deploy or move the targets on the ground during each scan. In previous studies, simultaneous measurement of illumination intensity with downwelling irradiance sensor on the ground or on-board platforms (i.e., UAV platform) that reflects illumination changes, was used to derive reflectance data [32]- [34], which does not provide a suitable solution(s).
The reflectance factor (used interchangeably hereinafter as reflectance) can be derived by the approach as follows [35], [36]: where λ is the wavelength, S is the spectral radiant intensity measured as the total amount of photons (DN) over ground objects (e.g., vegetation, soil) within the field of view (solid angle), D is dark current or noise measured with the lens cap on (unit is the same with S), and R is the spectral radiant intensity recorded by the sensor for white reflectance reference target, RC is the reflectance factor of the reference target. Thus, the reflectance factor of the sample image is also could be expressed as follows: where rfl img is the reflectance of experiment image (i.e., plants or soil), img smpλ is the sensor-recorded reflected spectral radiant intensity from the sample imagery, drk refλ is dark current/noise, which was measured periodically with the lens cap on. wht refλ and rfl whtλ are the sensor recorded upwelling spectral flux density (irradiance) of a Lambertian reference panel [37], and its reflectance factor, respectively. In our case, it is not practical to reliably deploy the Lambertian reference panel on the ground for each of the hyperspectral data cubes during the scanning. Thus, the relationship between the spectral flux density (flx dwnλ ) of a downwelling irradiance sensor and the upwelling spectral flux density from a Lambertian reference panel was characterized to extract the upwelling spectral flux density of the Lambertian reference panel from the downwelling irradiance sensor. To some extent, this approach essentially is a type of "continuous panel method," which provides reflectance factors by taking consecutive measurements of the sample target and Lambertian reference panel (converted from downwelling irradiance in our case) [38]- [40]. The conversion factor or function denoted as CF λ To characterize the CF λ , diurnal, as well as multiseasonal synchronized measurement of flx dwnλ using a downwelling irradiance sensor, and spectral flux density from a Lambertian reference panel using the hyperspectral scanner were conducted. A multistep Spectralon target with known spectral characteristics was used as a Lambertian reference panel, which was set up free of shadows or obstruction during the data collection. The two downwelling irradiance sensors (Ocean Optics, Largo, FL, USA) with a spectral range of 400-1000 nm (VNIR), and a spectral range of 900-2500 nm (SIWR) mounted on top of the field scanner measured stationary irradiance on board with 5 s intervals. CF is denoted as follows, and was modeled using an empirical line fitting approach: where f (·) denotes a function. Based on (2) and (3), the reflectance factor of each experimental imagery could be derived through the following equation:

C. BRDF Correction
The intrinsic anisotropy reflection characteristics of land surfaces including vegetation and soil, as well as various solar sensor viewing geometries within and between images, affect the spectral information of hyperspectral imagery, which often could be characterized by BRDF to generate normalized reflectance data [41], [42]. Studies have stated that radiance normalization and BRDF correction are essential for time-series applications and estimation of biophysical properties [34], [43]. We implemented a kernel-driven semi-empirical BRDF correction method rooted in the Ross Thick-Li Sparse (RTLS) model [44], [45] due to its proven advantage with numerous applications for images at different scales [46]. It contains three adjustable coefficients that are used for optimal BRDF fitting based on directional reflectance observations, and it characterizes surface reflection as a simple linear combination of three types of scattering: volumetric, geometric, and isotropic scattering [47]- [49]. The general formula is given as follows [45], [50]: where R is the surface directional reflectance, θ i is the incident zenith angle [generally referred to as the solar zenith angle (SZA)], θ v, is the reflected zenith angle [generally called the viewing zenith angle (VZA)], and ϕ r is the relative azimuth angle (Fig. 3). K vol is the volumetric scattering kernel (RossThick kernel) caused by multiple scattering within canopies, and K geo represents the geometric-optical scattering kernel associated with single scattering among canopies. f iso represents the contribution from isotropic scattering, f vol and f geo server as the coefficients/weights of the kernels at wavelength λ. A variety of kernel functions have been developed in previous studies. For geometric-optical scattering kernel (K geo ) [45], Li Sparse-Reciprocal kernel (LSRK) that assumes sparse canopy, Li-Dense Reciprocal kernel (LDRK) for a dense canopy [51], which often provides superior performance in the case of high solar and/or view zenith angles [48], [52], as well as sparse and dense transition Li Transit-Reciprocal kernel (LTRK) [53] were developed in previous studies. The kernels can be derived as follows [47], [51], [54]: The values for the parameters h/b and b/r (object shape and height) were determined after iterative testing [49].
The vector coefficients f iso , f vol, and f geo are extracted using multiple sample points based on linear function fitting and the least-square method [47]- [49]. As shown in (15), the anisotropy factor (ANIF) is defined as the normalization or ratio of reflectance data to fixed angle or nadir reflectance (user-defined specific geometrical angle, denoted as fixed), it is often used to describe the anisotropy characteristics of the surface target. Thus, the BRDF corrected reflectance rfl BRDF could be derived using (16) [49], [55], where rfl obs is the observed reflectance

D. Soil Detection and Removal
Background soil often affects the spectral characteristics of crops [56], which eventually reduces the quality and accuracy of remote-sensing-based crop monitoring and plant phenotyping [57]. Various pixel unmixing (separating soil and vegetation) methods have been attempted to reduce or minimize soil background effects on vegetation for coarse-resolution satellite remote sensing imagery. For instance, vegetation indices, MCARI [58], SAVI [59], and relative normalized difference vegetation index (RNDVI) [60], could differentiate dual peaks that corresponded to vegetation and soil reflectance.
These vegetation indices perform better over dark soils [61]. The combination of TCARI/OSAVI vegetative indices performed better for a closed canopy [62] while an MCARI/OSAVI combination performed better for an open canopy [63]. First-and second-order derivatives are also efficient in removing the effects of soil reflection spectra, although the second-order derivatives outperformed the former [61].
With the development of low-altitude aerial or proximal imaging systems, high-resolution remote sensing imagery becomes easily accessible, and detection and classification of soil and vegetation pixels deems necessary, which is of great importance in a variety of agricultural applications. Particularly, discrimination of soil and vegetation pixels also helps to estimate canopy cover, which is a significant agronomical component strongly related to crop growth, development, water use, and photosynthesis [64]. Estimation of the canopy cover area of crops may contribute to improving the efficiency of crop management practices and breeding programs [65].
A series of approaches have been developed to discriminate soil and vegetation pixels from high-resolution remote sensing imagery. The thresholding method is one of the commonly used approaches, such as automatic thresholding using the Otsu method [66]. In addition, RGB-based color indices-based threshold setting according to ratios of R/G, B/G [67], modified excessive green index (MExG), and excessive green index (ExG), have been employed to obtain soil-vegetation binary images [68]. A complex threshold using the combination of different indices has been applied as well [69]. Moreover, thresholding based on multispectral images has also been attempted [70]. For instance, NIR-based vegetation indices, NDVI, Green NDVI (GNDVI), and normalized green (NG), and associated performance over RGB imagery have been assessed as well [66]. The thresholding methods are relatively simple, but may not be sufficiently adaptive and robust for dynamic field environments and multitemporal cases [71], particularly images captured under various illumination conditions. Unsupervised or supervised classification or segmentationbased approaches have also been employed to discriminate soil and vegetation pixels, such as unsupervised method K -means clustering [72], as well as supervised classification method logistic function [73], support vector machine [74], and different segmentation methods [75], [76]. Classification and segmentation methods would provide accurate results, but are often time-consuming and computationally extensive. Additionally, supervised classification or segmentation methods often require a large number of training samples, which is also time-consuming and prone to human error [72]. Most importantly, none of these methods allow for fully automated processes for soil removal.
In this study, a fully automated and efficient method was proposed to discriminate soil and vegetation pixels from hyperspectral imagery that applies to both VNIR and SWIR spectral data. It is based on the physics of light interactions with soil and vegetation, i.e., photosynthetically active vegetation presents a reflection plateau at the NIR spectral region due to the cellular structure of leaves and strong absorption from chlorophyll in the blue and red wavelengths [77]. Thus, the reflectance of vegetation at the NIR region is higher than the VIS region, and the reflectance at the green region is higher than at the red and blue regions. The reflectance of soil (sunlit or shaded) demonstrates completely different trends from vegetation, gradually increasing from VIS to NIR regions (Fig. 4). Therefore, a vegetation and soil classification rule for VNIR data is derived as follows. Given a hyperspectral data cube I ∈ R m×n×b , where m × n is the image spatial dimension and b denotes the number of spectral bands. For an arbitrary pixel p i (where i is the index of the pixel and i ∈ {1, 2, . . . , m × n}), the proposed classification rule is expressed in the following: where the value 1 indicates photosynthetically active vegetation pixels whereas 0 indicates soil or dry vegetation pixels; B, G, R, and NIR represents blue, green, red, and NIR wavelengths, respectively. If the spectral profile of each pixel meets the criteria in (17), it will be classified as photosynthetically active vegetation, and if not, it will be regarded as soil or dry/dead vegetation. The hyperspectral SWIR sensor covers spectral ranges from 895 to 2501 nm. Due to the leaf cellular structure characteristics, along with the effect from vegetation biophysical quantities, vegetation (sunlit or shaded) often reflects strongly and exhibits a local peak at the spectral region around 1020-1120 nm (Region 1 ). In addition, vegetation also shows relatively lower reflection at the spectral region around 1160-1300 nm (Region 2 ) (Fig. 4), which mainly attributes to the absorption of leaf water [78], [79]. Thus, photosynthetically active vegetation demonstrates higher reflection at Region 1 than Region 2 , while soil (sunlit or shaded) often presents increasing reflection from Region 1 to Region 2 (Fig. 4). On the basis of spectral characteristics of vegetation and soil at these two spectral regions, if reflectance of a SWIR pixel at Region 1 is higher than Region 2 , the pixel is determined as a vegetation pixel. To this end, we proposed the following strategy/rule for vegetation/soil differentiation.
First, instead of using the average reflectance values of the two regions, we placed a stricter rule that selects the minimum reflectance value of Region 1 , denoted as Region1 Min and the maximum reflectance value of Region 2 , denoted as Region2 Max , which showed higher adaptivity to the noise of hyperspectral reflectance curve with superior performance. Additionally, a normalized approach was adopted as well, and a vegetation and soil classification rule, namely, VSDR SWIR , was proposed if V S D R SW I R > 0, vegetation pixels if V S D R SW I R < = 0, soil or dry vegetation pixels. Regardless of sunlit or shaded target, green band reflectance factor is always larger than those of blue and red and smaller than NIR reflectance factor for vegetation in VNIR region. For SWIR wavelengths, local maxima for region 1 (shown in gray columns) is always larger than local maxima of region 2 for vegetation and the opposite is true for soil.
For a given pixel v i in a SWIR hyperspectral data V ∈ R m×n×b , where m × n is the image spatial dimension and b denotes the number of spectral bands, if VSDI SWIR value is larger than 0, it will be classified as photosynthetically active vegetation, and if not, the pixel will be regarded as soil or dry/dead vegetation.
It is worth noting that, in order to reduce the impact from hyperspectral reflectance noise, which would potentially affect the performance or even lead to the failure of the above-proposed vegetation/soil pixel classification rules for VNIR and SWIR hyperspectral imagery, the Savitzky-Golay filter [80]- [82] was applied to the reflectance data to obtain smoother curves before the classification procedure.
The effectiveness of both VNIR and SWIR soil detection approaches was compared with commonly used soil masking techniques. In terms of VNIR images, we compared our proposed algorithm with an adaptive threshold-based [83] soil masking from NDVI, and a supervised support vector machine classifier trained from labeled pixels. The assessment was conducted for nine reflectance images (three images each from winter wheat, sorghum, and lettuce) collected from the field scanner. The image pixels were randomly labeled as vegetation or soil pixels and 70% of the data (n = 33 912) was used for training the model, whereas the remaining 30% (n = 14 535) was utilized for model evaluation. Since the number of features was very large (939 bands), a principal component analysis was applied to the training dataset with two components. Later, a support vector machine classifier was trained with tuned hyperparameters from the grid search algorithm. In the case of SWIR images, another support vector machine model was built using a similar procedure to that of VNIR images. However, the training set and testing set for SWIR images were 3051 and 1308, respectively, with 201 bands. Also, the principal component analysis for SWIR data was done with three components. Finally, after training the classifiers, the evaluation was done by calculating the overall accuracy, F1 score, and Kappa-coefficient for each method using the testing set from both VNIR and SWIR images.

E. Shadow Detection and Masking
Shadows prominently appear in high-resolution remote sensing imagery where objects (i.e., leaves or vegetative structures) block direct light produced by solar illumination [84]. In a remotely sensed reflectance imagery, shaded areas provide incomplete spectral information, lower intensities, and fuzzy boundaries, that lead to incorrect interpretation [85], [86]. The position of the camera box in the gantry hyperspectral system (Fig. 1) results in the camera box shadow being cast on all or a portion of plot imageries during specific times of the day, which varies by seasonal solar angles. Interpreting spectral values from imagery for plots with the shadow of field scanner may lead to inconsistent results, for instance, error in vegetation index calculation [87]. Therefore, we proposed a supervised classification system to detect the field scanner shadow from the hyperspectral images.
Multiple studies have utilized unsupervised classification, supervised classification, index-based methods, or physicalbased models to detect shadows. In this pipeline, we selected a supervised classification scheme because the goal was to detect shadows cast from the field scanner itself. The reflectance spectra from a shaded leaf and a sunlit leaf from the same hyperspectral imagery clearly showed distinct differences (Fig. 5). We developed a supervised classification scheme like one that is described in Section III-D for soil pixel identification; however, the training and testing samples for shadow detection were created separately (478 training samples and 206 testing samples). The pipeline included feature scaling, principal component analysis, and support vector machinebased classification.

F. Automated Image Quality Assessment
High spatial frequency components in an image can be attenuated due to focal blur or motion blur, which results in loss of details. To quantify the emergence of blur effect in hyperspectral VNIR and SWIR images, we propose to use a no-reference perceptual blur metric (NRPBM) where the value is ranging from 0 to 1 that, respectively, represent the best and the worst quality in terms of blur perception. The advantage of using the no-reference quality metric is that it neither requires prior knowledge of the original image nor assumption on the image content or the cause of the blur effect. This is particularly suitable for our case due to the fact that there is no standard reference available to compare in our database. The blur metric proposed by Crete et al. [88] was applied to the luminance component of an image, whereas in our case, it is applied to every spectral band of the input hyperspectral imagery. The algorithm measures the blur effect in an image by comparing the intensity variations between the original image with its blurred form. Fig. 6 shows a flowchart of the algorithm steps that refers to the following description.
Let I ∈ R m×n×b be a hyperspectral image with a spatial dimension of m × n pixels and a spectral dimension of b bands. The first step of the algorithm is to construct a blurred image which is achieved by filtering every spectral band, denoted as S j (the j th band for ∈ [1 b] ), with a low-pass filter in both horizontal and vertical directions, expressed by where F j H and F j V are the filtered j th band images in horizontal and vertical directions using filters h H and h V , respectively, given by where the superscript T denotes the transpose operation. Note that the filter h can be the other types such as Gaussian filters, averaging filters (e.g., 3 × 3, 5 × 5, and 7 × 7), or motion filters. Next, the intensity variations of the neighboring pixels are obtained by computing the absolute difference images as follows: where , y), and D F j V (x, y), and are the absolute difference images computed from the neighboring pixels in horizontal and vertical directions using the original input image and blurred images F j H and F j V . The variation of the neighboring pixels between the original and blurred images is then obtained by the following equations, which, respectively, indicate vertical and horizontal directions: To compare the variations from the original image, the sum of the values of D S j V , D S j H , V jd , and H jd is computed as follows: where subscript d refers to difference. Then, we normalize the above summation process to the range from 0 to 1 Finally, a no-reference perceptual blur estimation for j th band, denoted as B j , is obtained by During the quality evaluation process, the same process is applied to all bands of hyperspectral imagery separately, which are then checked for any abnormal values or for values of B j above a certain threshold such as 0.5.

G. Georeferencing
Geolocation of the hyperspectral data was accomplished through the conversion of scan gantry coordinates to UTM. The gantry x and y coordinates are referenced from a fixed position at the southeast corner of the field (0,0). As the camera box moves, the exact gantry coordinates representing its position over the plots based on the location in meters of the physical point at the southeast corner of the camera box are known and recorded as imagery metadata at the start of each camera scan. Additionally, the location of each sensor inside the camera box is known. The center pixel of the camera scan in the north-south direction is directly below the center of the sensor lens. Additional pixels in the north-south direction of each scan can then be calculated using pixel size for the specific sensor with half of the pixels to the north and half of the pixels to the south of the center pixel. The gantry is a line scanning system that moves in the east-west plane; hence, for each camera scan, the gantry location of every pixel in the east-west direction can be determined by adding or subtracting each pixel size to the camera scan's origin.
The specific latitude and longitude of the (0, 0) reference position were measured to be 33.074543 and −111.97479 with a known accuracy of 1 cm. Each image pixel is converted from gantry coordinates to UTM through conversion equations providing degrees per meter that were calculated specific to the latitude and longitude of the field location. The geolocation was then validated using targets with known gantry positions. The code-generated locations of the target positions in each hyperspectral scan that included a target were compared with the known positions of the targets. The code was determined to geolocate scans with an accuracy of +/−12 cm.

H. Automated Pipeline Architecture
All of the methods described above are applied to hyperspectral datasets automatically as part of a larger processing pipeline composed of several containerized components. When new datasets are collected from the gantry imaging system at the Maricopa Agricultural Center in Arizona, they are automatically transferred to the National Center for Supercomputing Applications (NCSA) at the University of Illinois using the Globus platform to move data and verify file integrity in the timestamp-based destination directory structure.
Due to the size of the datasets-a single file is normally 20 GB or larger and a full scan can total upward of 800 GB-it is critical to minimize the amount of repeat data transfer during processing. To that end, a network of virtual machines on NCSA's Nebula cluster each have direct access via a mounted drive to those data and execute the actual processing pipeline without having to download the files first. Data are uploaded to the Clowder research data management platform that acts as a web frontend for the TERRA-REF project, and Clowder in turn submits those datasets to the processing modules known as extractors automatically.
Each extractor is a Docker container running in the Nebula cluster that contains all scripts, software, and dependencies necessary to run the calibration algorithms and with the full storage system mounted and locally accessible. These containers can easily be duplicated and scaled as necessary to deal with incoming data flows. The primary limiting factor in scalability is the current algorithm's RAM requirement of 4× the raw input data size (e.g. 80 GB RAM for a 20 GB input file) which can be difficult to allocate with some of the larger datasets that are nearly 100 GB each. The extractor executes the calibration process and uploads metadata to Clowder summarizing the output. In the future, additional extractors can be added to the pipeline to further process those outputs and derive new insights.

I. Data Storage and Format
The continual capture, processing, and archiving of hyperspectral imagery present design challenges for the storage of each image, and ever increasing archive as a whole. Our design goals were threefold: 1) optimize datasets for automated harvesting of data to facilitate our own analyses via custom software, as well as web-originated queries through standard geoscience APIs and 2) annotate datasets with per-variable metadata that describes the meaning and units of fields, as well as storing provenance and processing metadata. A custom software was developed to convert imagery data to and then process it in a user-friendly, self-describing archival format that supports network transparent access by a wide variety of APIs. The workflow, which is based entirely on open-source software, annotates the datasets with comprehensive metadata information and propagated provenance information during each processing step.
First, the netCDF Operators (NCO) [89] converted the raw ENVI BIL format hyperspectral imagery to the netCDF4 enhanced data model format [90]. Images are stored as rectangular datacubes dimensioned (wavelength, y, x). This storage order optimizes the retrieval of single-wavelength hyperslabs through network data transfer protocols such as OpenDAP (DAP). The processing and archiving stream accesses these datasets via netCDF C, Python, (e.g., Web services), and OpenDAP bindings and protocols. Data are stored in the most compact datatype (e.g., uint16, float32) capable of representing the intrinsic precision of the data.
To further reduce storage (and data transmission) costs and bandwidth, the workflow may be optionally to create and store datasets using any lossy and lossless compression algorithms supported by netCDF4 [91]. Second, the datasets conform to the Climate and Forecast (CF) Metadata Conventions, version 1.8 [92]. Metadata are stored hierarchically as netCDF attributes in multiple top-level groups (e.g., gantry_system_fixed_metadata, sensor_fixed_metadata) to promote legibility. The data processing provenance is contained in the global history attribute. One novel feature stored in metadata is a pasteable text string that points the user to the exact location of the imagery using the Google Maps Engine.

A. Radiometric Calibration
The empirical relationship between hyperspectral VNIR and SWIR sensors recorded reflected spectral flux density from a Lambertian reference panel with known reflectance factor, and spectral flux (flx dwnλ ) measured synchronously by VNIR and SWIR downwelling irradiance sensors were established, respectively. Band-wise strong linear relationships characterized with high R 2 (the coefficient of determination) between hyperspectral sensor DNs and downwelling sensor irradiance values were observed for both VNIR (Fig. 7) and SWIR (Fig. 8) sensors (only a few typical bands were selected for demonstration here), and band-wise linear function that plays transformation role between DNs and irradiance values flx dwnλ was built for each band of both VNIR (Fig. 7) and SWIR (Fig. 8) sensors, and CF(flx dwnλ ) was computed and provided to the following equation that was mentioned in the methodology section, to derive reflectance values of VNIR and SWIR hyperspectral imagery: It is worth noting that, the data points used to compute the CF λ were collected diurnally, as well as at multiple seasons, to account for low, median, and high solar geometry. Additionally, both VNIR and SWIR downwelling irradiance sensors have higher spectral resolution, thus spectral resampling was applied to VNIR and SWIR downwelling spectral flux data to match with a spectral resolution of VNIR and SWIR hyperspectral sensors, respectively. Moreover, to collect one hyperspectral cube data often takes about 5 min in our case, however, synchronized downwelling irradiance sensors have a much higher frequency which takes measurements at the second level, thus, the corresponding irradiance records within the time frame of each hyperspectral data cube collection (i.e., 5 min) were averaged and used in the calibration procedures. Last but not least, compared to SWIR data, VNIR data presented a stronger relationship between spectral flux (also known as DNs) from Lambertian reference panel and downwelling spectral flux (flx dwnλ ) with higher R 2 , as well as more convergence pattern of data points (Figs. 7 and 8), which is likely due to sensitivity/response of the sensor array to the reflected flux at different wavelengths, as well as the noise level at different wavelengths caused by the SNR of  VNIR and SWIR sensors [93]. In future work, more complex and nonlinear, as well as machine learning-based conversion methods should be evaluated.
Based on the above-mentioned procedure, VNIR and SWIR calibration models were developed and applied for VNIR and SWIR hyperspectral imagery over different crops. Fig. 9 demonstrates the profiles of both raw values (also known as DNs) and reflectance of VNIR imagery from sorghum, durum wheat, and soil, respectively. As shown in Fig. 9, the reflectance of both sorghum and wheat leaves displayed typical vegetation spectral profile characteristics, which present strong reflection at the NIR region, followed by the green spectral region, and strong absorption at the red and blue spectra regions [94], [95]. Additionally, the reflectance values at different wavelengths also shown typical ranges, for instance, reflectance values at the NIR region varied around 0.6-0.7 ranges. Moreover, the reflectance of soil exhibits a typical soil spectral profile pattern as well, which often increases from VIS to NIR spectra regions [96]. It is worth noting that, regardless of whether a pixel was sorghum, durum wheat, or soil, reflectance values at the region of wavelengths higher than 900 nm display a relatively noisy pattern, which is due to the lower SNR of hyperspectral VNIR sensors at longer wavelengths [93].
For SWIR hyperspectral data, wavebands with high noise often due to water absorption (atmospheric windows) [97], [98] at the regions of 1306-1437, 1790-1992, and 2400-2500 nm were dropped. The reflectance of sorghum, durum wheat, and soil displayed a typical pattern of vegetation at the SWIR spectral region (Fig. 10). Strong reflectance Fig. 9. VNIR spectral profiles of sorghum, durum wheat, and soil before and after radiometric calibration. Reflectance factor values at different wavelengths fall within typical data ranges observed for vegetation and soil. One may expect that vegetation reflectance at NIR region to be 0.55-0.6 and our results varied around 0.6-0.7 ranges. This might be attributed to uncertainties resulting from exposure time and scanning speed for this sample image, BRDF effect, or both. Reflectance of soil exhibits a typical soil spectral profile pattern as well, which often increases from VIS to NIR spectra regions.
around the wavelength of 1100 nm due to biophysical quantity and yield, and slight absorption around the region of 1180 nm due to leaf water content, as well as two "peak shape" region around 1500-1700 and 2000-2400 nm often caused by reflection/absorption effect due to plant lignin, cellulose, and water content [79], [99], [100]. The reflectance of soil also showed a typical pattern of soil spectral profile at the SWIR region [101].
To further assess the calibration results through a quantitative manner, the VNIR and SWIR calibration models were applied to hyperspectral images that include a Lambertian Spectralon reference panel. The multistep Spectralon reference panel has known reflectance values of 99%, 50%, 25%, and 12%. It is worth mentioning that those testing images are independent of the reference images used for calibration model building. As shown in Fig. 11, the red, blue, green, and pink lines represent averaged spectral profiles from the regions of the panel with 99%, 50%, 25%, and 12% reflectance factors, and the black continuous lines represent the ground truth reflectance factors of each region stated by the reference panel vender. Overall, the reflectance profile ([ Fig. 11(b) and (c)] of each region of the multistep Spectralon panel from the radiometrically calibrated image [ Fig. 11(a)] presents a good matching with the ground truth reflectance, indicating good performance of the calibration models. For VNIR data, spectral profiles at the spectral region with wavelength larger than 950 nm demonstrated overestimation and lower convergence pattern with the ground truth reflectance profiles. This is likely due to the lower SNR and higher noise characteristics of hyperspectral sensors at longer wavelengths in the NIR region [93]. For SWIR data, spectral profile from 12%, 25%, and 50% panel regions showed a higher convergence trend with ground truth values, while spectral profile from 99% panel region presented a lower convergence pattern. For further analysis, root-mean-squared error (RMSE) between calibration results and ground truth reflectance values was computed (Fig. 12). Compared to the VNIR calibration model, the SWIR calibration model yielded lower RMSE for all the four 99%, 50%, 25%, and 12% panel regions, indicating a superior radiometric calibration performance. It is worth noting that, the RMSE at 12% and 25% panel regions are lower than 50% and 99% panel regions, implying the SWIR calibration model performs better for targets with lower reflection characteristics. The VNIR calibration model obtained the lowest RMSE value for the 50% panel region, implying it produces the best results for the targets that have median reflection patterns.

B. BRDF Correction
The efficacy of BRDF correction was evaluated by analyzing VNIR and SWIR hyperspectral image reflectance data  collected at different solar and sensor viewing conditions before and after BRDF correction. As shown in Figs. 13 and 14, for both VNIR and SWIR data, before and after BRDF correction, the reflectance of the same spot [average value of the yellow rectangle region in Fig. 13(a)] from five diurnally collected images were compared, the results showed that the BRDF corrected reflectance profiles tend to be more convergent, particularly the NIR spectra region are closer to typical vegetation spectral characteristics in terms of reflectance values. It is worth noting that, the convergent pattern of SWIR data after BRDF correction is not as obvious as VNIR data [Figs. 13 and 14]. This was also quantitatively evidenced by band-wise standard deviation (STD) for data samples before and after BRDF correction (Fig. 15). The STD of both VNIR and SWIR data samples are reduced after BRDF correction (Fig. 15), which indicate that the BRDF correction was able to  VNIR spectral patterns of the same target location at five different times from morning to afternoon before and after BRDF correction. (a) True-color visualization of the test image, (b) spectral profiles before BRDF correction, and (c) spectral profiles after BRDF correction. Vegetation spectra of the same pixel from diurnal images become more aligned to one another, demonstrating the improvement of BRDF correction. SWIR spectral patterns of the same target location at five different times from morning to afternoon before and after BRDF correction. (a) False-color visualization of test image, (b) spectral profiles before BRDF correction, and (c) spectral profiles after BRDF correction. It is evident from the results that SWIR reflectance becomes consistently similar after the BRDF correction. normalize the differences potentially caused by various solar incident and azimuth angles from morning to the afternoon to some extent, and provide more temporally comparable reflectance values.  The reflectance values of pixels located at different viewingangle regions of the hyperspectral imagery (i.e., center and two sides of an imagery) were compared before and after applying BRDF correction models. As shown in Fig. 16, after BRDF correction, the VNIR spectral profiles from the left, middle, and right positions of the imagery (which represent different viewing angles), demonstrated slightly convergent patterns, which was also shown by the decreased STD at most of the spectral regions [ Fig. 18(a)]. Fig. 17 presents the differences in spectral profiles from pixels at the left, middle, and right positions of SWIR imagery before and after applying the BRDF correction model, which showed a slight convergent trend after the correction as evidenced by the slight decrease in STD [ Fig. 18(b)].
The convergent pattern displayed in our case after BRDF correction to reflectance values with various sensor viewing angles is consistent with previous studies [48], [102], indicating the BRDF correction is able to produce reflectance  data that are comparable at various viewing angles. It is worth noting that, compared to the diurnal data samples, after BRDF correction, the change of STD of reflectance values derived from pixels at different viewing angles is less apparent (Figs. 15 and 18). This might be due to the fact that the change in viewing geometry at different viewing angles is much smaller compared to the change of solar incidence and azimuth geometry at a diurnal time scale.

C. Soil Detection and Removal
The proposed rule-based soil masking method for VNIR images yielded the best performance (94.17% accuracy) followed by the SVM-based (93.92% accuracy) supervised method and the NDVI-based (91.87% accuracy) threshold method (Table II). The performance metrics were calculated for around 15 000 independent pixels sampled from multiple VNIR images. Fig. 19 shows the results after applying three different masks to three crops, i.e., durum wheat, sorghum, and lettuce. The average NDVI for each image before and after removing shadow is also shown. The NDVI is an indicator of green healthy vegetation, where higher numbers represent more vegetation [103]. The average NDVI increased for all Fig. 19. VNIR images after application of different soil masking methods, i.e., rule, NDVI, and SVM-based models, for durum wheat, sorghum, and lettuce plants. The average NDVI of each image is provided which represents the overall vegetation. For each method and plant, the average NDVI increases after applying the soil mask. The NDVI and SVM-based methods failed to accurately detect leaf edges as vegetation pixels compared to the unsupervised rule-based method. The rule-based method developed for SWIR images also showed promising results in detecting soil pixels. The supervised SVM-based method slightly outperformed (92.58% accuracy) the rule-based method (91.21% accuracy), which was conducted for around 5000 independently sampled pixels from SWIR images (Table II). The effect of removing soil pixels from SWIR images is visualized in Fig. 20 along with the corresponding average normalized difference water index (NDWI). The NDWI was calculated using the reflectance from the 1080 and 1700 nm wavelengths as discussed in Ghulam et al. [104], where higher values of NDWI represent more canopy representation. The average NDWI increased for each method and crop (e.g., durum wheat, sorghum, and lettuce) after removing the soil. However, the rule-based method proposed for SWIR images inaccurately classified soil pixels as vegetation pixels which resulted in Fig. 20. SWIR images after applying rule-based and SVM-based soil masking methods for durum wheat, sorghum, and lettuce. The false-color composite was created with reflectance from 1596.86, 1197.44, and 1032.97 nm as RGB. The average NDWI for each image is provided, which represents the average canopy water content before and after soil removal. The unsupervised rule-based soil removal method performed closely with supervised SVM-based method as the average NDWI value is similar. slightly lower accuracy, and NDWI value compared to the SVM-based method.
The rule-based method proposed for VNIR outperformed other commonly used soil removal techniques, whereas the SWIR rule-based method reached similar performance to a supervised method. The simplicity, flexibility, and applicability of the unsupervised rule-based methods are promising since they do not require any training sample generation or model tuning which significantly reduces computation complexity.

D. Shadow Detection and Masking
The shadow detection models for both VNIR and SWIR images were specifically trained to identify shaded regions cast by the gantry hyperspectral system [ Fig. 1(b)]. The SVM-trained classifier for VNIR and SWIR yielded 96.74% and 97.55% accuracy, respectively, which was carried out using independently selected sample pixels. Fig. 21 shows sample visualizations of shadow masks for different crops. Both shadow detection models could distinguish between shaded regions, sunlit leaves, and sunlit soil for different crops. The automated pipeline included VNIR and SWIR images with corresponding shadow masks.

E. Automated Image Quality Assessments
The hyperspectral image quality in terms of sharpness or other artifacts (e.g., motion blur) related to the blur can be measured by blur metrics such as NRPBM mentioned in Section III-F. To examine the effectiveness of NRPBM, four different typical images in our dataset were selected based on crop growth stages and shaded/sunlit conditions. Fig. 22(a)-(d) show RGB color composite images for visual demonstration and Fig. 22(e) shows the corresponding band-wise NRPBM value of each hyperspectral sample images. It can be observed that there are significant differences in terms of loss of details (blurriness) of leaves, between the Image-4 (most clear and high quality) and the other three images (i.e., Images-1, 2, 3) with varying degree of blurriness. This is echoed with the band-wise NRPBM values where the NRPBM score of Image-4 [ Fig. 22(d)] is the lowest (which Fig. 21. Demonstration of SVM-based gantry shadow detection model with (a)-(d) VNIR and (e)-(h) SWIR images. Two separate models trained for VNIR and SWIR accurately detect shaded pixels generated by the field scanner. Fig. 22. Examples of RGB color bands in the hyperspectral dataset used for image quality measurement (images were stretched using 5% linear method for better visibility) and corresponding image quality measurement using NRPBM (lower value indicates better quality). Among the sample images, (c) Image-3 has the lowest image quality because of shadow and (d) Image-4 has the highest quality (lowest values) compared to other images as evidenced by sharp leaf edges and clear ground objects. means better) across all spectral bands compared to that of the other images. Although Image-4 exhibits sharper edges on the leaves of sunlit areas, its shaded regions on the leaves result in lower visibility than Image-1 [ Fig. 22(a)] and Image-2 [ Fig. 22(b)] therefore affect the NRPBM score as evidenced in higher NRPBM values. Comparing, their visual quality is similar and NRPBM scores are comparable, especially for the wavelengths shorter than 720 nm. The comparative graph in Fig. 22(e) also indicates that the spectral band quality, in terms of NRPBM, of all four images tends to be gradually improved after the wavelength cross approximately 784 nm, while continuously decreases between 400 and 528 nm. This may be due to the sensor characteristics in those spectral bands related to mechanical sensor design and light dispersion, as well as physical hardware calibration of the sensor. It is worth mentioning that blur effect and shaded/sunlit conditions are not only seen in RGB bands as shown in Fig. 22, but also in other corresponding spectral bands. As expected, Image-3 [ Fig. 22(c)] showed the poorest image quality (highest NRPBM values over the entire wavelengths) due to a large portion of the imagery is covered by shadow.

V. CONCLUSION
Here we presented a fully automated, data-driven approach to the calibration of super high-resolution hyperspectral big data. Following a rigorous radiometric calibration, BRDF correction and reflectance normalization, soil and shadow masking, and image quality assessments were carried out automatically. The pipeline also derives selected plant traits, including spectral vegetation indices, canopy cover, and stress quantification. The results demonstrated that the proposed workflows are valid and effective for handling hyperspectral big data, thus providing an additional path to accelerate plant breeding and high-throughput crop phenotyping for improved yield and biomass traits through advanced imaging technologies.
1) Comparison of spectral profiles from calibrated hyperspectral imagery of a multistep (12%, 25%, 50%, and 99%) Spectralon panel with vendor-provided ground truth reflectance values demonstrated that RMSE between the ground truth and calibrated reflectance factor is about 0.05-0.075 for VNIR data and 0.025 for SWIR data. Compared to the VNIR calibration model, the SWIR calibration model yielded lower RMSE for all the four 12%, 25%, 50%, and 99% panel regions. The RMSE at 12% and 25% panel regions are lower than 50% and 99% panel regions for the SWIR region, implying that the SWIR calibration model performs better for targets with lower reflection characteristics. The VNIR calibration model provided the lowest RMSE value for the 50% panel region, implying it produces the best results for the targets that have median reflection patterns (e.g., vegetation). This might also be attributed to the fact that radiometric calibration models were developed using a 50% reference panel which is closer to vegetation reflectance values. 2) BRDF correction was carried out to normalize the data over a diurnal cycle and across sensor view angle differences from nadir to the edge of each scan line. The BRDF correction significantly improved both VNIR and SWIR reflectance factor data demonstrated by a closer alignment of reflectance profile of a ground target over diurnal cycle. Band-wise STD for the five diurnal samples collected from morning to afternoon and over different pixels of the same target over the range of viewing angles was reduced up to 50% after BRDF correction, indicating improvement. 3) A fast, rule-based method for soil and shadow removal was proposed to automate the process. This method was compared against commonly used machine learning algorithms, including support vector machine as well as vegetation index-based thresholding. Our proposed rule-based approach outperformed conventional methods with up to 95% accuracy observed for the VNIR and SWIR data. Established on the foundations of remote sensing physics, this fast and effective method holds great promise for the automation of big data processing as it does not require any training samples. 4) In order to automatically detect issues like blur effect, mechanical failure of the scanning mirror, or changes in the SNR from environmental disturbances, we proposed an NRPBM method for hyperspectral image quality assessment. This pipeline produces a quality metric between 0 and 1, respectively, representing the best and the worst quality. This allowed the data-driven calibration system to automatically detect issues and trigger timely warnings for events that required physical investigation.