A co-aligned HySpex (Norsk Elektro Optikk, Norway) system was utilized to capture hyperspectral imagery during each field season. This system was made up of a VNIR-1800 camera with 186 spectral bands between 400 nm and 1,000 nm and a spectral resolution of 3.26 nm, as well as a SWIR-384 camera with 288 spectral bands between 953 nm and 2,518 nm and a resolution of 5.45 nm. The flight was at an altitude of approximately 365 m above the ground, giving a spatial resolution of 0.25 m in the VNIR-1800 system and 0.5 m in the SWIR-348 system. Imagery collection dates (as indicated by crops’ days after emergence (DAE)) in 2020 and 2021 were chosen to ensure they covered the critical potato growth stages in each field season. They included: early tuber bulking (53 DAE in 2020 and 51 DAE in 2021), mid tuber bulking (61 DAE in 2020 and 59 DAE in 2021), and late tuber bulking (81 DAE in 2020 and 85 DAE in 2021). Imagery collection dates were coordinated with field sampling dates as closely as possible, although weather factors (primarily precipitation and cloud coverage) limited the ability to align flight dates perfectly. For predictions of in-season traits, ground truthing data was paired with imaging data only when their collection dates were no more than 3 days apart. Pairs of collection dates for aerial imaging and in-season field measurement in this study are shown in Table 2. Data from Russet Burbank and Soraya individually and combined were used for in-season predictions. End-of-season traits were predicted using at-harvest measurements paired with spectral data collected at early, mid- and late tuber bulking separately in each of the 2 years. In addition, models for end-of-season traits were developed using spectral data collected from all three growth stages combined across the 2 years. This paper reported at-harvest predictions for Russet Burbank because only this variety was grown in both years.

HySpex images were processed using a consistent and open-source workflow following Liu et al. (2021), https://github.com/EnSpec/hypro). Processing steps included radiometric calibration using manufacturer-provided coefficients, correction for smile effects (shifts in wavelength position across spectrometer focal array (Schläpfer and Richter, 2010), orthorectification to <1 m precision using GPS-IMU data, atmospheric correction (Richter and Schläpfer, 2010), and correction of bi-directional reflectance distribution effects due to variations in sun-sensor-target geometry across flightlines (Queally et al., 2022).

Spectral smoothing and data extraction for statistical analyses Savitzky–Golay filter (window size: 7) was used to reduce noise in surface reflectance (Schafer, 2011). Then, the average reflectance of all vegetation pixels within each subplot was extracted. Non-vegetation pixels (e.g., shadows and soil) were excluded using NDVI<0.5. Vector normalization was performed to reduce the effects of illumination conditions on reflectance (Asner et al., 2015; Feilhauer et al., 2010). Reflectance values extracted from the hyperspectral imagery were averaged over each plot for modeling. Spectral data from atmospheric water absorption features at 1,320–1,470 nm and 1780–2000 nm were removed from analysis, as were noisy data from 2,350–2,500 nm. Plot mean reflectance values were paired with the target predicted trait measured from each plot and were then used to build partial least squares regression (PLSR) models. All steps were performed using the open-source Geospatial Data Abstraction Library (GDAL, version 3.0.4), as well as Python packages Pillow (6.2.0), Numpy (1.17.2) and Scipy (1.3.1).

All PLSR modeling was done using the HyTraits package in Python. PLSR analyses (Figure 1) consisted of an outer loop, where the original dataset was split into a calibration (70%) and testing (30%) dataset, and an inner loop, where the calibration dataset was further split into a training (70%) and validation (30%) dataset, which was used for hyperparameter testing. Models were built using different values for the maximum number of components in the PLS process, with the model that performed the best on the validation data being saved. The outer loop was nested within the inner loop, with each inner loop process repeated 100 times, and a total of 200 outer loop processes. This pipeline yielded 200 final models that were saved from each model’s respective inner loop process. Final predictions were created by averaging the output of these 200 final models. Models were evaluated according to the testing R ² and Normalized Root Mean Square Error (NRMSE). The standardized regression coefficients were used to evaluate the contribution of each wavelength to the estimation of potato traits (Liu et al., 2021).

Cross-year validation was conducted where the models were calibrated using data from 1 year and tested on data from the other year. Only data from Russet Burbank that was cultivated in both years was used for the cross-year validation. Cross-variety validation was performed using data from one variety for calibration and those from the other variety for testing. In addition to models using the entire spectrum (400 nm–2,500 nm), models were built using only the visible-near infrared (VNIR) spectrum (400 nm–1,300 nm) or the shortwave infrared (SWIR) spectrum (1,400 nm–2,500 nm) (Liu et al., 2021).

All data management was performed using the Python (version 3.8.8) packages NumPy (version 1.20.1), pandas (version 1.2.4), and SciKit Learn (version 0.24.1) and the R (R Core Team, 2019) package ggplot2 (version 3.4.2).

For analyses using both years, model performance across all five in-season traits of interest and by variety exhibited R ² > 0.5 and NRMSE ranging 0.15–0.27 (Figure 2). Compared to the control treatment, the three higher N rates (Rate one–3) were consistently associated with higher N status indicators, and higher total N uptake in both aboveground whole canopies and belowground tubers. However, no general trends could be summarized about N effects on the in-season immediate tuber yield. In most cases, the models created using the Russet Burbank data alone slightly outperformed those created using the Soraya data alone or models created from both varieties.

Models built from the leaf data had better performance (total leaf N, R ² = 0.66–0.70) than those built from the petiole data (petiole nitrate-N, R ² = 0.54–0.58). Models predicting whole canopy N uptake showed moderate accuracy, whereas those predicting tuber N uptake and in-season immediate yield were strong with R ² values ranging from 0.52 to 0.71.

Prediction of traits measured at harvest are shown in Figure 3. The control treatment consistently produced the lowest yield, smallest tuber size, and lowest fertilizer N recovery % at harvest. No trend of N effects on specific gravity were observed. Marketable yields in 2021 were higher than those in 2020 across all treatments, likely due to different weather conditions during the two growing seasons. In general, models built from imaging data collected from the three different growth stages had comparable results, with the best results from early- and mid-tuber bulking stages. In all cases, models predicting marketable yield and specific gravity all performed well (R ² > 0.6), whereas models predicting fertilizer N recovery % exhibited lower evaluation statistics (R ² < 0.4).

Models to predict at-harvest traits using spectral data combined from all three growth stages (named multi-stage models) for Russet Burbank typically did not improve in comparison to those built using only one date of imagery (named single-date models), except for marketable yield (data not shown). No general trends emerged when comparing multi-stage models with the single-date models developed using imagery from the mid- or late tuber bulking stages.

In general, models calibrated with Soraya and tested with Russet Burbank generated similar or better prediction results (similar or higher R ²) for both aboveground and belowground traits than those calibrated with Russet Burbank and tested with Soraya (Table 3). The belowground tuber N uptake could be well validated both ways with R ² close to or higher than 0.6. Comparatively, predictions of the immediate tuber yield could not be well validated across the two varieties.

For the cross-year validation, the multi-stage model predicting marketable yield showed better performances with higher R ² and lower NRMSE than the single-date models for both ways of validation (Table 4). Specific gravity could not be well validated in either case. Fertilizer N recovery % showed moderate validation results across different growth stages or over the whole season.

Full spectrum standardized regression coefficient values (“coefficients” in short) illustrate the contribution of different wavelengths or wavelength regions to models predicting in-season aboveground and belowground traits (Figure 4). For interpretation, these are used to identify the plant spectral characteristics that contribute to the prediction capacity, in particular since some traits (e.g., marketable yield, specific gravity) are not directly sensed in the imagery, but rather are inferred from foliar characteristics. In both the VNIR and SWIR regions, coefficient weights typically followed consistent trends among the different trait prediction models. The coefficients of features from the regions (530–630 nm and 700–750 nm) mostly associated with chlorophyll had relatively large weights. Coefficients from the red-edge position (680 nm–780 nm) were weighted heavily as well. Although coefficients from the near-infrared region were not consistently weighted heavily, some portions of this region had high weights, particularly from 900 to 1,200 nm (the region that is related to plant structure).

Features had overall lower coefficients in the SWIR region than in the VNIR region. And features in the SWIR region showed more difference between the two varieties, which is likely related to their genetic differences. In general, features at the beginning (from 1,500 nm to 1,600 nm) and the end (from 2,300 nm to 2,450 nm) of the SWIR region had larger weights. Coefficients from the SWIR region were overall the highest in models built from the Soraya data alone and were the lowest in those built from the Russet Burbank data alone. Over the whole spectrum, although following the same trends, coefficients at specific wavelengths of Russet Burbank are different from those of Soraya, indicating that those two varieties have different spectral signatures under the same N treatments.

The standardized regression coefficient values for models using the images collected at one of the three growth stages across the 2 years for Russet Burbank to predict at-harvest traits follow the same trends as those found in the in-season belowground traits, with the highest weights found in the 900 nm–1,200 nm range (Figure 5). No consistent trends were noted across different traits for models built from imagery taken at different growth stages (i.e., early vs. mid-vs. late tuber bulking). For models predicting marketable yield and fertilizer N recovery % by tubers, the coefficients of features from the chlorophyll-associated wavelengths were weighted lower at the late tuber bulking stage of the season than in the early season, whereas the opposite was true for models predicting specific gravity. Models built using imagery from the mid-tuber bulking stage had larger weights on features in the 2,250 nm–2,450 nm range than those built using imagery from the early or late tuber bulking phases.

Predictions of each in-season trait were performed using three different spectral datasets: 1) the entire set of spectral data; 2) the visible-near infrared (VNIR) region ranging from 420 nm–1,300 nm; and 3) the shortwave infrared (SWIR) region ranging from 1,400 nm–2,350 nm. The best spectral region (with the highest R ² and lowest NRMSE) for each trait of the two varieties individually and combined is reported in Table 5. The use of different spectral regions generally resulted in comparable R ² and NRMSE values for each trait, but predictive models built from the two varieties individually or combined had noticeably different performances under different spectral regions. For Russet Burbank, the models using the SWIR spectral data consistently underperformed those using the VNIR spectral data or the full spectrum. In comparison, when creating predictive models with only the Soraya data or with data of both varieties, no single spectral region consistently performed the best. Predictions of plant N status using petioles and leaves, and the whole canopy N uptake generally favored the VNIR region.

Similar to in-season predictions, evaluation statistics were often similar between models predicting the at-harvest traits using different spectral regions (Table 6). Models predicting fertilizer N recovery mostly showed better performances using only the VNIR region, while those predicting tuber marketable yield were less consistently impacted by spectral regions, for both single-date and multi-stage models. Notably, the SWIR region at the late tuber bulking stage produced the best prediction results for all the three at-harvest traits. For multi-stage models, using the full spectrum always produced the best prediction results.

For the predicted indicators to assess crop N status (petiole nitrate-N and leaf total N), higher N rates consistently resulted in higher crop N status indicators. Predictive models built from petiole nitrate-N showed moderate R ² and NRMSE (Figure 2). Although petioles are the primary plant parts that transport inorganic N in the form of nitrate within the plants, and nitrate-N levels are the most common indicators to guide potato farmers’ in-season N fertilization decisions (Colquhoun, et al., 2020; Zebarth and Rosen, 2007), agronomic studies often see high levels of variability in petiole nitrate-N levels under varied N treatments, as they can be influenced by many factors such as time of tissue collection during the day, cloud cover, and temperature (Liu, et al., 2021; Wang, et al., 2022). In contrast, total N content in the leaves does not exhibit this high variability, and, as well, foliage (in contrast to petioles) is directly observed by the sensor. With a higher number of field measurements, the variability in petiole nitrate-N might be reduced, possibly resulting in stronger predictions. While additional research is needed to further refine the understanding of the relationship between plant tissue types and the spectral data collected from them, the results of this study suggest that whole leaf total N (directly observed) are better than petiole nitrate-N for monitoring and predicting potato N status, and should be used instead of petiole nitrate-N to inform in-season N fertilization if remote sensing technologies are applied (Liu, et al., 2021).

Predictive models of in-season N uptake by both aboveground canopies and belowground tubers of each and both varieties had moderately good performance with R ² values between 0.43 and 0.64 (Figure 2). Tuber initiation leads to a sharp increase of N absorbed and utilized by the plants, and tuber bulking signals a transition for potato plants to rapidly partition N and energy from the aboveground canopy to the belowground tubers (Alva et al., 2002; Lauer, 1985; Millard et al., 1989; Zebarth and Rosen, 2007). Typically, 70%–85% of the N in the plants is present in the tubers, although as much as 40% of plant N can be present in potato vines (Zebarth, et al., 2004). Our results suggest that hyperspectral remote sensing of the foliage could be moderately effective in estimating potato N uptake both aboveground and belowground. Notably, predictions of tuber N uptake were slightly better (with higher R ²) than those of the whole canopy N uptake (Figure 2). This capability likely results from a tight linkage between aboveground N (detected in imagery from leaves) and N availability belowground to tubers (which is not directly imaged). This correspondence can be seen in the standardized PLSR coefficients for whole canopy total N (Figure 4, third row) and fertilizer N recovery % (Figure 5, third row), which are comparable and show higher coefficients in similar wavelength ranges.

In addition, in-season immediate tuber yield predictions had favorable evaluation statistics. Photosynthetic activity is one of the major factors affecting tuber growth, and spectral reflectance data can reflect photosynthetic activity well (Gitelson and Merzlyak, 1996; Gitelson and Merzlyak, 1998; Patane and Vibhute, 2014). Besides, plant vigor is generally highly correlated with crop yield and is reflected well in hyperspectral images due to high reflectance of starch, cellulose, lignin, and nonstructural carbohydrates in the near-infrared and shortwave infrared regions (Curran, 1989; Liu et al., 2021). These factors may contribute to the good real-time prediction for the immediate tuber productivity using the hyperspectral signatures of the foliage. Note that standardized PLSR coefficients for intermediate yield and nitrogen-related traits track closely (Figure 4), further supporting the interpretation that the ability to predict yield depends on the ability to characterize factors related to photosynthesis, such as plant N.

Models predicting end-of-season marketable yield performed well (R ² = 0.77–0.78). Final marketable yield could be well predicted as early as the beginning stage of tuber bulking (R ² = 0.77, NRMSE = 0.2), and the good prediction results could be found throughout the mid- and late tuber bulking periods, showing that higher N rates could produce higher marketable yields (Figure 3). Accurate prediction of marketable yield from in-season imagery over the course of the growing season is greatly valuable to potato farmers as it could inform economic and management decision-making for producing sellable agricultural products, support sustainable practices that underscore efficient N fertilizer use, and reduce impact of crop production on the environment and the groundwater quality (Chang, et al., 2023; Luczaj and Masarik, 2015).

Models predicting at-harvest specific gravity showed good evaluation statistics as well (Figure 3). As a major indicator of potato crop quality, predictions of end-of-season specific gravity complement the good marketable yield predictions, allowing for farmers to further estimate the final value of their crops before the end of the season so they can adjust their harvest time and practices to ensure the optimal conditions for maximum crop profitability. The models derived from late-season imagery did not perform as well as those derived from early and mid-season imagery (R ² = 0.61 compared to 0.74 and 0.68). Tuber starch accumulation, a quality parameter that tuber specific gravity indicates, tends to level off during the later part of the growing season than in the early or middle of the season, when tuber bulking occurs at an exponential rate (Zebarth and Rosen, 2007). This may result in a relatively weaker late-season prediction for specific gravity of the tubers.

Fertilizer N recovery is the ratio of total amount of N that leaves a cropping system through the harvested tubers to total seasonal applied fertilizer N. The N that was not recovered by tubers at harvest will be lost to the groundwater via nitrate leaching, to the atmosphere via denitrification, and to terrestrial and aquatic systems through wet or dry deposition, all of which could cause environmental concerns (Zebarth and Rosen, 2007). Compared to the better performance of models predicting the in-season tuber N uptake (Figure 2), evaluation statistics for models predicting fertilizer N recovery at the end of the season were relatively weak (R ² < 0.4, Figure 3). This could suggest that unlike for other potato traits, hyperspectral remote sensing might not have the capacity to predict a trait that needed to be calculated compared to a trait that could be directly measured.

For cross-variety testing, our study showed that models calibrated with Soraya, a Dutch variety bred to be a highly efficient N user, could be moderately tested on Russet Burbank, the No.1 variety developed and grown in the U.S. for over 100 years that is values for its frozen fry processing quality but requires high nutrient and water inputs. Our findings suggest that the new elite potato varieties possess similar in-season traits as the traditional ones, but otherwise they harbor some unique characteristics such as the high N use efficiency and low need of fertilizer inputs, which could make them more sustainable to grow. There is ongoing research to replace Russet Burbank in the U.S., and varieties such as Soraya from other germplasms could provide good genetic resources to facilitate that effort. This indicates that the development of new varieties bred for different characteristics will be best predicted from hyperspectral imagery using models that explicitly include those varieties in the training data.

For cross-year testing using Russet Burbank, we did not achieve satisfying results to predict the at-harvest traits from the single-date models, but the multi-stage models to predict marketable tuber yield produced relatively strong results. This finding demonstrates that there are systematic differences in plant growth caused by different temperature and precipitation conditions in each season (Supplemental Figure 1), or in the hyperspectral images collected at different growth stages of each year. These differences may not be accounted for in the single-stage or single-season models to predict specific gravity and fertilizer recovery %. However for marketable tuber yield, if we combined images over all three growth stages (early, mid-, and late tuber bulking) that covered the growing condition range throughout the season, it could generate moderate testing results between the 2 years. In conjunction with the cross-variety results, this result points to the need to develop training data sets that cover the reasonable variation that would be expected across varieties, growth stages and years.

The PLSR coefficients were similar to those observed in Liu et al. (2021) for petiole nitrate-N, whole leaf total N, and yield at harvest. Positions of the most important wavelengths in predicting the in-season traits were consistently found in the visible and red-edge region (530–630 nm and 680–780 nm) that is highly associated with chlorophyll absorption and canopy greenness, the near-infrared region (900–1,200 nm) that is directly related with the canopy structure and overall vegetation health (Brewer, et al., 2022; Croft, et al., 2020), and at the beginning (1,500–1,600 nm) and the end of the shortwave infrared region (2,300–2,450 nm) that are related to plant starch, cellulose and protein content (Fathi, 2022). Although the positions of the most importance wavelengths stayed within the same ranges over the spectrum between the two varieties (Figure 4), specific wavelengths were typically differed between Russet Burbank and Soraya, noting that all the in-season traits measured in this study were related to plant N status and N uptake. This finding again suggests that although the foliar hyperspectral signatures generally showed similar relationships with plants’ biophysical and biochemical characteristics, some wavelengths were still variety-specific, caused by the unique N use efficiency and partitioning between aboveground and underground N of that variety.

When comparing models built from the VNIR region alone, the SWIR region alone, or the full spectrum, the best spectral regions that generated the highest R ² and lowest NRMSE were trait- and variety-specific, again pointing to differences in plant spectral and nutrient usage characteristics resulting from genetic differences. Models predicting the aboveground traits (N status and whole canopy N uptake, Table 5) consistently performed better when the VNIR region was used. Specifically for whole leaf total N prediction, when using the full spectrum to build models, the contribution from SWIR wavelengths seemed to be suppressed so that using the full spectrum and only using the VNIR generated the same model evaluation statistics, possibly due to very small variations in the SWIR reflectance among different N treatments (Liu, et al., 2021). For belowground traits (tuber N uptake and immediate tuber yield, Table 5), the best regions are highly variable between different varieties/combination of varieties. The VNIR region contains wavelengths where plant spectral signatures are impacted by chlorophyll activity, likely leading to the observed effects of VNIR on the photosynthetic aboveground tissues than on the belowground traits. In comparison, the ubiquitous preference of the SWIR region by Soraya models to predict the belowground traits (Table 5) may be influenced by its high N-use efficiency and special characteristics to partition N between the canopy and the tubers (D. Maum, personal communication, 10 January 2020), as the SWIR wavelengths have important features related to plant nutrient uptake and overall plant physiology (Gold, et al., 2020).

Table 6 clearly shows that during the early and mid-tuber bulking phases, the best regions to predict traits at harvest mostly included the VNIR part of the spectrum, related to vegetation greenness and vigor, but they switched to the SWIR region during the late tuber bulking phase, more related to foliar nutrition and foliar chemical content. Results from our PLSR models using the late tuber bulking images indicated the importance of the SWIR region to predict potato traits at harvest. It has been reported that SWIR wavelengths are related to crop dry matter content. During the latter part of the growing season, potato dry matter content accumulation plateaus and is maximized, while the aboveground canopy is senescing and wilting (Zebarth and Rosen, 2007), and therefore the utility of VNIR data that’s related to chlorophyll and canopy structure would be substantially reduced, and the SWIR-only model became the best option for predicting potato traits using late-growth-stage imagery. As well, information from the SWIR region may be more strongly related to nitrogen absorption features associated with proteins (Berger, et al., 2020), despite of confounding effects from chlorophyll and canopy vigor. When making predictions on data that include aerial images across different growth stages of the crop, using the whole spectrum including both the VNIR and SWIR regions may result in more robust predictive results (Table 6), as the relationships between plants’ spectral reflectance and their biophysical and biochemical properties changed over the course of the growing season.

Optimally managing nitrogen for potato production is extremely challenging. The amount of N the crop needs and the available amount in the soil change throughout the growing season. Ideally, N should be added as it is needed, matching nutrient supply with crop demand. Unfortunately, growers cannot currently make real-time nitrogen budgets and calculate the ideal amount of N to add at any certain time. Instead, they use guidelines developed through research and past experiences for fertilization scheduling. The most commonly adopted method in the U.S. is petiole nitrate-N testing, and the measurements are typically made every 5–10 days. Petiole nitrate-N suffers from high variability, high labor costs and long turnaround testing time (typically 3 days). The present study shows that by using narrow-band hyperspectral imagery collected at different potato growth stages, it is possible to develop relatively robust PLSR models with higher R ² and low NRMSE values for the most popular potato variety (Russet Burbank) in the U.S. and a new variety (Soraya) from Europe. The prediction results could help farmers to gain a real-time picture of their crop N status and yield forecasting, so that better N fertilization decisions can be made, unnecessary N use can be reduced, and the associated groundwater contamination may eventually be mitigated. To further refine this approach, additional spectral and field measurement data gathered from different locations and years, potato varieties, and production practices should be used for training and testing. It will create more robust models, capture further nuance behind variations in the relationship between spectra and crop biophysical and biochemical traits, and generate more finely tuned precision resource management tools.