Predicting Phosphorus and Potato Yield Using Active and Passive Sensors

: Applications of remote sensing are important in improving potato production through the broader adoption of precision agriculture. This technology could be useful in decreasing the potential contamination of soil and water due to the over-fertilization of agriculture crops. The objective of this study was to assess the utility of active sensors (Crop Circle ™ , Holland Scientific, Inc., Lincoln, NE, USA and GreenSeeker ™ , Trimble Navigation Limited, Sunnyvale, CA, USA) and passive sensors (multispectral imaging with Unmanned Arial Vehicles (UAVs)) to predict total potato yield and phosphorus (P) uptake. The experimental design was a randomized complete block with four replications and six P treatments, ranging from 0 to 280 kg P ha − 1 , as triple superphosphate (46% P 2 O 5 ). Vegetation indices (VIs) and plant pigment levels were calculated at various time points during the potato growth cycle, correlated with total potato yields and P uptake by the stepwise ﬁtting of multiple linear regression models. Data generated by Crop Circle ™ and GreenSeeker ™ had a low predictive value of potato yields, especially early in the season. Crop Circle ™ performed better than GreenSeeker ™ in predicting plant P uptake. In contrast, the passive sensor data provided good estimates of total yields early in the season but had a poor correlation with P uptake. The combined use of active and passive sensors presents an opportunity for better P management in potatoes.


Introduction
Potato (Solanum tuberosum L.) is among the most important crops grown in the United States of America (USA) [1]. Maine is ranked ninth among the states in area dedicated to potato production; however, yields are much lower than in Western and Midwestern USA [2]. Phosphorus (P) management is a crucial component of potato production because it impacts early root development and crop yield [3][4][5]. P plays a vital role in the plant's energy metabolism, such as the synthesis of adenosine triphosphate (ATP) and adenosine diphosphate (ADP) [4,6]. P is also an essential component of phospholipids, phosphoproteins, nucleotides, and coenzymes [7]. It regulates enzymatic metabolism the most widely used vegetation index for determining plant health based on the red wavelength at 660 nm and NIR (700-1200 nm) [36].
Studies on the application of remote sensing techniques to detect P deficiencies are scarce. Leaf reflectance (LR) in the blue band of the visible region (VR) of the electromagnetic spectrum is an excellent measure to estimate the P content in plants [15]. Although P stress had a higher reflectance in the spectrum's yellow and green portions, the red edge did not show a regular shift in chlorophyll absorption at the near-infrared (NIR) region [5,37].
It is crucial to consider the diversity of P spectral signatures of plants to perceive the variability in the signature [38]. The red-edge band (680-760 nm) can be used to predict leaf phosphorus concentration [39] accurately. The NIR region is more sensitive to the plant P status than the visible region [10]. While P stress decreases the LAI, there is no consistent pattern among the different P treatments [40]. The electromagnetic spectrum regions affected by phosphorus were in the range of wavelengths between 400 and 1300 nm [15]. Reflectance at the blue band (445 nm) and NIR (730-900 nm) regions were most suitable to predict P at the early growth stages [11].
The development of new technologies and an improved understanding of the cause-and-effect relationships regarding spatial variability could lead to better management of the applied nutrients, making such measures more economically and environmentally sound [15,35,41]. In this context, the objectives of this study were to compare active and passive sensor performance in detecting P deficiency, total yield, and P uptake in potato throughout the growing season.

Site Description
Twelve sites were set up on potato farms in Aroostook County, Maine, USA, in 2018 and 2019 (Tables 1 and 2). The potato cultivars, soil types, and preceding crops are shown in Table 1. The average temperature was between 12 and 20°C and the precipitation was between 50 and 80 mm. The experimental design at all sites was a randomized complete block with four replications and six P treatments: 0, 56, 112, 168, 224, and 280 kg P ha −1 applied as triple superphosphate (46% P 2 O 5 ). In 2018, each experimental unit (plot) was 9.1 m in length and 3.6 m in width, and the rows were spaced at least 0.9 m apart. Nitrogen (N) and potassium (K) were applied equally for all P treatments, with 202 kg N ha −1 as urea (46% N) and 224 kg K ha −1 as KCl (61% K 2 O). In 2019, the same P application rates and resources were applied, and each experimental unit (plot) was 6.1 m in length and 1.8 m in width; the rows were spaced at least 0.9 m apart. Nitrogen was applied equally for all P treatments, with 224 kg N ha¯1 as ammonium nitrate (34% N), K was applied equally for all P treatments as in 2018. The fertilizers were mixed and applied by the replacement method at planting. Pesticides, such as chlorothalonil, and herbicides, such as Weedar64 (Nufarm) 2,4-D, sourced from Alsip, Illinois, USA, were applied by the farmers as recommended by the University of Maine, Orono, Maine, USA [42].  Modified Morgan Extraction (MME); a P extracted by MME (P-MME); b Mehlich 3 (M3); Mehlich 3 soil test P (P-M3); c K extracted by M3 (K-M3); d Iron extracted by MME (Fe-MME); e Iron extracted by M3 (Fe-M3); f Aluminum extracted by M3 (Al-M3); g Organic Matter (OM); h Cation exchangeable capacity extracted by MME (CEC − MME).

Soil Sampling and Analysis
Before planting, soil samples were collected from the twelve sites (Table 1) using a 2-cm diameter stainless probe. There were five subsamples collected at a depth of 0 to 25 cm for pH, organic matter, P, N-NO 3 − , N-NH 4 + , and K analysis ( Table 2). The soil samples were air-dried and ground to pass through 2-mm sieves. Soil pH was estimated using a 1:1 soil: deionized water solution and a pH meter (pH robot is a labfit AS3000; company HQ located in Perth, Australia), with an appropriate electrode (s) [43]; the organic matter was determined using loss of weight on ignition [44]. The phosphorus (P) was analyzed by Mehlich 3 extraction, followed by analyzing the filtrate via the Inductively Coupled Plasma (ICP) (Spectro Genesis, Kleve, Germany) method [45]. The P and K levels were determined by Modified Morgan Extraction (MME) and estimated within the linear range on the ICP. The N-NO 3 − was measured by using an automated ion analyzer and the recommended ratio of 1:10 (5 g soil; 50 mL extractant (0.01 M CaSO4, 2 M KCl)) [46]. The N-NH 4 + was measured via KCl extraction [47].
The University of Maine Soil Laboratory, Orono, Maine, USA, analyzed the soil samples.

Plant Sampling and Analysis
One petiole and leaflet (whole leaves) was collected, at the tuber initiation stage (60 days after the planting date), from each potato plant, thus, four plants were collected from each sub-plot for the phosphorus analysis. The petioles and leaflets were taken from the fourth and fifth fully expanded leaves from the plant top [48], dried in an oven with air circulation at 65°C for 96 h, and ground. Subsequently, the samples were turned into ash at 550°C for 5 h in a muffle furnace. The obtained ash was dissolved in 50% HCl on a hot plate and analyzed using ICP [49]. The University of Maine Soil Laboratory analyzed the plant samples. Phosphorus uptake (kg ha −1 ) was calculated by multiplying the dry matter of the petioles and leaflets by the percentage of phosphorus content (%) as follows: Phosphorus uptake (kg ha −1 ) = Dry matter of petioles and leaflets (kg ha −1 )× Phosphorus content (%) 100 (1)

Potato Harvest
At the harvest time (after 130 to 140 days from planting; Table 1), the tubers were dug using a two-row digger and then handpicked manually from a 5.5 m 2 area within the central rows of each plot. Individual plots were harvested separately to remove possible border effects. Fresh tubers were graded and weighed to determine the total potato tuber yield (Mg ha −1 ).

Active Sensors
We used the sensors GreenSeeker™ (GS; Trimble Navigation Limited, Sunnyvale, CA, USA) and Holland Scientific Crop Circle™ ACS 430 (CC; Holland Scientific, Inc., Lincoln, NE, USA). The GS sensor was used to measure incident light and reflected light from plants at 660 ± 15 nm (red) and 770 ± 15 nm (NIR) to obtain the NDVI and IRVI [50]. Sensor readings were downloaded by connecting the sensor to a laptop computer. The GreenSeeker™ (GS) emits two wavelengths of reflectance light: 660 ± 15 nm (red band) and 770 ± 15 nm (NIR band). The GS active sensor light would be emitted from diodes in alternate bursts, thus, the red and NIR bands pulse for 1 ms (1 millisecond) at 40,000 Hz. Each burst from the diodes would obtain up to 40 pulses before pausing for the other diode to emit its radiation (plus 40 pulses). The lighted area is 60 cm wide by 1 cm long, with the extending dimension basically positioned vertically to the direction of the GS sensor in the field movement.
The CC sensor was simultaneously used to measure the crop/soil reflectance at 670 nm, 730 nm, and 780 nm. The data were collected using a Holland Scientific GeoSCOUT GLS-400 (Holland Scientific, Inc., Lincoln, NE, USA) datalogger and downloaded to a laptop computer. The CC yielded several vegetation indices, such as NDRE (normalize differences red-edge), NDVI, and the leaf area index (LAI), obtaining as the sensor output chlorophyll red edge (CHLRE) ( Table 3). The CC-430 active sensor incorporates three fixed wavebands. The field of view (FOV) was an oval of~30 • by~14 • range. The readings would be collected at one measurement height (60 cm) above the potato canopy at 10 Hz (10 readings/s) when moving at a constant speed in each experimental unit. The sensor path was parallel to the potato rows with the beam of reflected light positioned vertically to the potato rows.
Each sensor recorded about 45-60 readings during one passing through the plot and collected readings from the two middle rows. The reading values were averaged and organized using the in-house Excel macro option for Visual Basic (Microsoft Excel 2020, Version 16.43, NC, USA). For June to August, thirteen active sensing dates took place in each of the twelve experimental sites: 25 June, 1 July, 9 July, 12 July, 18 July, 22 July, 25 July, 1 August, 5 August, 13 August, 16 August, 20 August, and 23 August. The sensing dates were excluded from the results when there was no significant model fit to the data.

Passive Sensors
Unmanned Aerial Vehicle Image Acquisition and Processing In 2018, two DJI Phantom 4 UAVs (Shenzhen, Guangdong, China) ( Figure 1a) were used for image acquisition. One had a NIR camera/sensor and was used to obtain the NIR images ( Figure 1b). The other had a regular camera/sensor used to collect the visible-color images of the red, green, and blue bands (Figure 1c). The gimbal controlled the UAV movement (pitch and roll) throughout the flight [51]. In 2019, A DJI Inspire 2 UAV (Shenzhen, Guangdong, China) with a portable Altum multispectral sensor (MicaSense, Seattle, WA, USA) was used to collect the multispectral data ( Figure 1d). Figure 1e shows how the Altum sensor integrates a radiometric thermal camera with five high-resolution narrow bands and produces advanced thermal multispectral and high-resolution imagery in one flight for the advanced analytic spectral bands: blue (475 nm center, 20 nm bandwidth), green (560 nm center, 20 nm bandwidth), red (668 nm center, 10 nm bandwidth), red-edge (717 nm center, 10 nm bandwidth), NIR (840 nm center, 40 nm bandwidth), and thermal (LWIR thermal infrared 8000 to 14,000 nm, radiometrically calibrated). The altitude for the UAV was 75-80 m; the flight duration was based on the field size, which was between 8 and 12 min; the resolution was 3.2 cm pixels −1 ; the overlap was 90/90; the image count was 86 to 112 images, depending on the distance of the field; and maximum speed was 2.7 m/s. The flight was performed between 10:00 a.m. and 2:30 p.m. local time in the study area. A certification of authorization was obtained from the United States Federal Aviation Administration (US-FAA) in Lincoln, Nebraska, USA, for the DJI Phantom 4 and DJI Inspire 2 flights. In 2018, images were taken by a Sony 4K camera (Shenzhen, Guangdong, China) with a 12 megapixel 1/2.3" complementary metal-oxide-semiconductor (CMOS) sensor (4000 × 3000 pixels). Both camera settings were maximum wide-angle identical to a 28 focal length (35 mm format equivalent), with auto exposure bracketing, auto speed, autofocus, auto white balance, f/2.8 aperture, and focus at infinity ∞.
Acquisition points were generated in the field and used as a waypoint for the flight path. In each flight operation, a steady ambient lighting condition was attempted (sunny days). In all flight campaigns, Map Pilot for DJI (Man Made Easy, Shenzhen, Guangdong, China) was used. Maps and waypoints were generated and saved for the next flights. For the Altum sensor, a calibrated reflectance panel was used to calibrate images.
Image calibration was implemented before each flight using an automatic panel detection mode that automatically recognizes and obtains the calibration panel image. The live view tab of the passive sensor configuration page displays a live image and detects the reflectance panel and quick response (QR) mode after selecting the QR button using the sensor camera. The drone-sensor was held at least 1 m above the calibrated reflectance panel to obtain the image from centering the panel and QR code in the field of view. Once the sensor detected the panel, it sent a unique sound and flashed the blue LED light, indicating that the image was captured. To ensure that the calibrated image was taken, the calibrated panel should be placed flat on the ground and far from any objects that could impact the light and illuminate the captured image. The operator (drone pilot) stands in front of the calibrated panel, and the sunshade is directed at the operator's back. Shadows should be eliminated when the panel is on the ground to prevent invalidated reflectance compensation readings.
The experiment's images were uploaded via Agisoft Metashape software (Agisoft LLC, St. Petersburg, Russia) once for all bands. Photos were uploaded, including the reflectance calibration images. The calibrated images were selected and prompted to upload the calibration from the comma separated values (CSV) file provided by MicaSense (MicaSense, Seattle, WA, USA); this CSV file has all the calibration values corresponding to the MicaSence radiometric panel, obtained from the MicaSense website on https://micasense.com/prv/. Therefore, in the Agisoft Metashape processing software, the reflectance option and the sun sensor in the calibrated reflectance dialog were selected to perform the calibration based on panel data and image meta information.
Drones had the Downwelling Light Sensor (DLS), which provides reliable and accurate radiometric data. In this study, the drones' DLS provided a sun sensor that measures the solar irradiance at the drone's top. It is an important element to minimize the change in the environmental light conditions. The DLS normalized the data to the standard conditions.
An orthomosaic was created from images in 2018 using the Pix4DMapper software (Pix4D mapper, Lausanne, Switzerland) [52]. The Agisoft Metashape professional software 1.6.3 [53] was used to process the images from 2019. Rectifications to the camera's view distortions were performed, and the images were stitched together to create a geo-referenced orthomosaic (GeoTIFF). The GeoTIFF image for each flight was uploaded in ArcGIS (Version 10.3; Esri Headquarters, Redlands, CA, USA) [54]. In ArcGIS, polygons were generated to represent each experimental plot ( Figure 2). All processing was performed with an HP laptop (Spring, TX, USA), with a 4-core Intel 10TH Gen i7-1065G7 CPU @ 1.30 GHz (max 1.5 GHz), 16 GB memory, 1TB SSD + 32GB Optane, and NVIDIA GeForce RTX 2070 GPU(Palo Alto, CA, USA). The experimental field at the early leaf stage: (a) black boxes indicate experimental plots, and white boxes within the black boxes indicate the buffer area for each plot and for two middle rows, manually; (b) boxes indicate NDVI for all pixels and plots; (c) manually, non-plant pixels were removed from the image using reclassification; (d) non-plant pixels were excluded by multiplying (c) by (b) to obtain the NDVI image just for plant pixels; (e) NDVI values for each plot (average of two middle rows) in the field, in this photo, zonal statistics (spatial analyst) in Arc GIS is shown in this photo to differentiate the experiment plots, automatically. In our study, we used zonal statistics as a table to obtain the means values. The same procedure was done for all passive vegetation indices and plant pigments.
For June to August, eight passive sensing dates took place in each of the twelve experimental sites: 25 June, 1 July, 9 July, 16 July, 22 July, 29 July, 14 August, and 23 August. Some flight operation dates were excluded from the results when no model provided a significant fit to the data.

Visible Bands and NIR Vegetation Indices
The digital numbers (DN) converted reflectance was calculated by converting the *.TIF files to float files (*.FLT) (Figure 3). The VI equations in Table 3 were used as input to the raster calculator to calculate the VIs. To acquire the average for each plot (two middle rows), we used the command "Zonal statistics as a table". These processes were repeated for all the VIs and each orthomosaic photo (GeoTIFF) by employing an iterator.

Statistical Analysis
Because the sensing dates were the same for both years, the data were pooled for each date. As interpreted from the handheld sensors and the drones, both the active and passive sensors showed some skewing degree, indicating the potential benefits of data diagnostics and transformation, such as a log transformation. These analyses were performed on both transformed and untransformed data. All the data were analyzed using RStudio software (Version 1.3.1073; Boston, MA, USA).

Criteria for Evaluation of a Subset of Predictor Variables
Both untransformed and transformed models were evaluated by adjusted R 2 (R 2 adj ), Akaike's Information Criterion (AIC), Corrected Akaike's Information Criterion (AICc), Bayesian Information Criterion (BIC), and Prediction Sum of Squares (PRESS), as follows: where R 2 adj represents the adjusted R 2 ; SSE is the error sum of squares; y i is the ith observation; n is the number of observations; y is the mean of the n observations; SST the total sum of squares; y is the predicted value of y; and p is the number of predictors. R 2 adj was implemented in R using the R sq.adj package.
Akaike's Information Criterion AIC AIC = (n log (SSE/n) + 2p + other terms (5) where n is the number of observations; SSE is the error sum of squares; and p is the number of predictors. Other terms are independent of the error sum of squares and p. The smaller AIC values indicate the better model fit.
The AIC package was implemented in R software as AIC < -sapply (1:m, function (x) round (extractAIC (om [[x]], k = 2) [2], 2)). The other terms were independent of SEE and P and the same for every model.

Corrected Akaike's Information Criterion AICc
AICc is developed to overcome the model complexity when the AIC is not strong enough. In the multiple linear regression model, AICc is obtained as shown in Equation (6): where n is the number of observations and p is the number of predictors.

Bayesian Information Criterion BIC
The BIC is often used to select variables in multiple regression problems. In this study, BIC was used to determine the best model (where regression models were a variable subset). The following equation can obtain the BIC: where K = p + 2, the number of parameters estimated in the regression model; and MLE is the maximum likelihood estimation used to estimate the probability distribution parameters by maximizing a likelihood function.

Prediction Sum of Squares (PRESS)
The PRESS summarizes the multiple linear regression model's fit to a sample of observations using the fitted regression function to obtain the predicted value. PRESS is shown in Equation (8): where y i is the predicted value for the ith subject; and y i − y i is the prediction error for the ith subject.

Generalized Linear Model (GLM)
For the GLM, we first tested the overall correlation between the total potato yield and the vegetation indices for active and passive sensors. In the untransformed model, we selected the best predictors (p < 0.05) and excluded the non-significant predictors to increase the adjusted R 2 and minimize the subsection variables; the multiple linear regression approach was used with the transformed model. Backward and forward stepwise regression was used with all vegetation indices for the active and passive sensors to guide which model might have the minimum AIC and BIC in predicting total potato yield and phosphorus uptake. Finally, models were established based on the maximized adjusted R 2 to minimize all possible selections. These vegetation indices were used to predict the total potato yield and phosphorus uptake based on the GLM, as shown in Equation (9): where Y i indicates the response variable; X in represents the predictor variables; β 0 , β 1 , and β n are the parameters of the model; and ε i is the error term.

Crop Circle™ Total Potato Yield and Phosphorus Uptake Models
Based on the adjusted R 2 , p value, AIC, AICc, BIC, and PRESS values (Table 4), the most affected total potato yield variables were NDRE, NDVI, and CHLRE.
The models for the first sensing date, 25 June, and the last significant sensing date, 1 August, showed the highest adjusted R 2 values (0.36 and 0.31), respectively. The last sensing date (1 August) model was improved after transforming (R 2 adj = 0.32). The lowest R 2 adj values (0.10) were obtained for the fourth, sixth, and seventh sensing dates, 12 July, 22 July, and 25 July, respectively. At the first sensing date (25 June), the multiple regression model applied to the test data set had the best R 2 adj (0.36) with the smallest AICc and BIC values. Several models were improved following the transformation. All the final regression models were statistically significant (all p values < 0.0001), except for 22 July's untransformed model ( Table 4).
The P uptake models had the best fit on the third sensing date, based on the highest R 2 adj value (0.60) and the lowest validation criteria ( Table 5). The fourth and fifth sensing dates showed a weak correlation between the response variable (P uptake) and predictors (NDVI, NDRE, LAI, and CHLRE). The NDVI was the only vegetation index that impacted P uptake on the fourth sensing date (12 July), with an R 2 adj of 0.32. There were no significant responses for the last six sensing dates; therefore, the equations were excluded from Table 5.

Relationships between the Actual and Predicted Variables
The strongest relationship between the actual and predicted total potato yields was observed for both the transformed and untransformed models on 25 June (Figure 4a,b) and on 1 August (Figure 4e,f). In contrast, the weakest relationship was found for the untransformed model on July 22 (Figure 4c). However, transforming the data increased the variation explained by the model (Figure 4d).
The strongest relationship between the actual and the predicted P uptake was detected on 9 July for both the transformed and untransformed data (Figure 5a,b). In contrast, the weakest relationship was foreseen on 25 July (Figure 5d). The transformation had little effect on the model fit on any of the sensing dates ( Figure 5).

Total Potato Yield and Phosphorus Uptake Models
There was little correlation between the observed and predicted values until August 5 ( Table 6). The transformations increased the coefficients of determination values on 5 August, 16 August, and 20 August (Table 6).
On the contrary, the models best predicted phosphorous uptake early in the season (Table 7).

Relationships between the Actual and Predicted Variables
The most robust relationship between the actual and predicted total potato yield was obtained by the model using the transformed data on August 5 (Figure 6b). The weakest relationship with total potato yield was found for an untransformed model on 16 August (Figure 6c). For the P uptake, the coefficient of determination was comparable to those obtained for potato yields. The most robust relationship was obtained for the untransformed model on 1 July (Figure 7a). Transformations did not have a major input on model performance.

Total Potato Yield and Phosphorus Uptake Models
Models constructed using the data from the passive sensors provided a relatively good prediction of the total potato yield throughout the growing season, starting on the earliest sensing date of 25 June, with coefficients of determination fluctuating around 0.50 (Table 8). The vegetation indices that most frequently predicted total potato yield were NDVI, CHLGR, and ANTHO. Transforming the data did not result in better-fitted models (Table 8). Total potato yield was negatively correlated with ANTHO on the first drone flight date (June 25). P uptake models explained very little variation observed in our experiments (Table 9). Once again, transforming the data did not produce models with meaningfully higher coefficients of determination. The outcomes were consistent across all sensing dates (Table 9).

Relationships between Actual and Predicted Variables
Models built using the data obtained by the passive sensors showed stronger and more consistent correlations between the actual and predicted total potato yields than the models based on the data from active sensors. The strongest relationship between the actual and predicted total potato yield was obtained by the transformed model for the first drone flight on 25 June, with an adjusted R 2 of 0.63 (Figure 8b). For other models, the coefficients of determination fluctuated between 0.50 and 0.56. On the other hand, the phosphorous uptake data showed little correlation between the actual and predicted values (Figure 9).

Discussion
Understanding the relationships and correlations between yield and phosphorous application requires assessing different sites with different vegetation indices and plant pigments [11]. The NIR region plays a vital role in predicting phosphorus stress due to the internal leaf structure; under P stress, the number of small leaf cells is increased compared to non-stress conditions [11,58].
Indices obtained by the Crop Circle™ active sensor were weakly correlated with total potato yield, most likely because the Crop Circle™ active sensor can be more sensitive to nitrogen stress than to P stress [12,59]. The highest amount of variation explained by the model was only 36% (R 2 adj = 0.36), achieved at the first sensing date (25 June). After that, the fitted models' explanatory capabilities dropped (Table 6). In contrast, 60% of the variation in P uptake (R 2 adj = 0.60) could be explained by the model on the third sensing date (9 July). However, subsequently, sensor performance deteriorated ( Table 7). This could be attributed to potato tissue phosphorus that could be decreased at the end of the season growth stages due to the P translocation from the shoot system to tubers regardless of the soil P concentration, which could exhibit a vast P stress [15,60].
Correlations between the actual and predicted variables revealed the poor predictive abilities of Crop Circle™ regarding the total potato yield, as well as being only a moderate predictive model for P uptake. This could be attributed to the values of NDVI, which were lower than those at NDRE due to the red band, which was not useful to predict phosphorus stress and the overlapping between plant pigment reflectance [15].
For the GreenSeeker™ active sensor output, the adjusted R 2 values were weak to moderate for NDVI and IRVI at all dates. This might be because the GreenSeeker active sensor operates in two wavelengths centered at the red and NIR with no red-edge. The red-edge band can predict P in crops [61]. The GreenSeeker™ active sensor readings showed a weak relationship with potato yield and P uptake when used individually. Similar to Crop Circle™, it was a useful sensor for nitrogen rather than phosphorus [12,62].
In this study, the vegetation indices derived from passive sensors showed a better correlation with total potato yield than the active sensors, which might be due to the extraction and analysis of the images based on removing the soil and other non-plant pixels before calculating the vegetation indices using the passive sensors [63,64]. At the first flight, the ANTHO and NDVI were negatively correlated with total potato yield, with a moderate adjusted R 2 (0.50), and positively correlated during the third drone flight, with an adjusted R 2 of 0.56. ANTHO was estimated from the red band and green band's reflectance ratio to reduce the overlap between the ANTHO and chlorophyll [17]. However, there are two possible explanations for this result. Firstly, the leaf surface reflectance in the visible region is larger than that in the infrared region. Secondly, the ANTHO content is high at the early vegetative stage, characterized by low photosynthetic activity [13]. At the early and late leaf stages, the green band might have a lower reflectance (higher absorption) than the middle stage of the growing season [14,21]. P uptake was not successfully determined by the models created based on the passive sensor data. During the flights, the adjusted R 2 values ranged from 0.06 to 0.17. Even though the transformed methods reinforced the equations, the models were still considered poor predictive models. However, early potato stages, such as vegetative growth (30 days after planting) and tuber initiation (50-60 days after planting), are important stages for the plants to accumulate sufficient soil P, which can be reflected in potato growth [6,60]. At the end of the potato growing season, the cell number, shape, size, and water content decrease, leading to increased absorption at the green band and increased blue and red reflectance. Furthermore, plant stress at the end of the potato growing season could be confounded with plant senescence [18]. Additionally, there is no huge P stress at the late potato growth stage due to the high soil P, thus, P stress would be rarely detected by the sensors [15].

Conclusions
This study assessed the employment of active sensing (Crop Circle™ and GreenSeeker™) and passive sensing (multispectral imaging with UAVs) to predict total potato yield and P uptake. The vegetation indices were calculated from the data collected by the active and passive sensors to generate the models. The vegetation indices developed using Crop Circle™ showed weak correlations with total potato yield but were good in providing an early season estimate of P uptake. Similar results were observed using GreenSeeker™. GreenSeeker™ is also moderately effective in predicting P uptake early in the season. The passive sensor data extracted from visible and multispectral imagery provide good estimates of total potato yield at the early potato season and can be used to predict P stress though it predicted P uptake only late in the season. The remotely sensed imagery showed a significant correlation between total potato yield and the NDVI, CHLGR, and ANTHO indices. However, it was largely useless in predicting P uptake. Future investigations should assess these results on soils with different P concentrations, under different environmental conditions, and using different active and passive sensors.
Author Contributions: A.J., L.K.S., A.A., C.W. and S.K.B. designed, collected and analyzed data, and wrote the manuscript. A.Z. and A.B. collected soil, and plant data; L.K.S. and A.A. acquired the funding. All authors have read and agreed to the published version of the manuscript.