Nitrogen Estimation for Wheat Using UAV-Based and Satellite Multispectral Imagery, Topographic Metrics, Leaf Area Index, Plant Height, Soil Moisture, and Machine Learning Methods

To improve productivity, reduce production costs, and minimize the environmental impacts of agriculture, the advancement of nitrogen (N) fertilizer management methods is needed. The objective of this study is to compare the use of Unmanned Aerial Vehicle (UAV) multispectral imagery and PlanetScope satellite imagery, together with plant height, leaf area index (LAI), soil moisture, and field topographic metrics to predict the canopy nitrogen weight (g/m2) of wheat fields in southwestern Ontario, Canada. Random Forests (RF) and support vector regression (SVR) models, applied to either UAV imagery or satellite imagery, were evaluated for canopy nitrogen weight prediction. The top-performing UAV imagery-based validation model used SVR with seven selected variables (plant height, LAI, four VIs, and the NIR band) with an R2 of 0.80 and an RMSE of 2.62 g/m2. The best satellite imagery-based validation model was RF, which used 17 variables including plant height, LAI, the four PlanetScope bands, and 11 VIs, resulting in an R2 of 0.92 and an RMSE of 1.75 g/m2. The model information can be used to improve field nitrogen predictions for the effective management of N fertilizer.


Introduction
Precision agriculture (PA) is a management technique that selectively applies crop farming resources such as fertilizer, water, pesticides, and herbicides based on the plant needs within a field [1][2][3]. Nitrogen is an essential macronutrient to plants as a major constituent of organic material including enzymic processes, chlorophyll, and oxidationreduction reactions; levels of nitrogen in plant tissue can indicate yield potential and crop health [4]. However, nitrogen is one of the most expensive nutrients to supply, and studies found that nitrogen recovery efficiency by annual crops was, on average, less than 50% of the amount of fertilizer applied [5,6]. Excessive fertilizer can leach from the soil and contaminate waterways, disrupting local ecosystems and causing denitrification that results in greenhouse gas emissions [7]. Nutrients that have been added beyond the critical level of maximum growth can continue to accumulate in the plant tissue without any further yield increase [4]. Commonly in grain crops such as wheat, excessive nitrogen can cause plant stems to grow tall to the point of lodging-the stems bend over, making it difficult to harvest and increasing the chances of grain moisture and disease, and often reducing yield significantly [8]. Usually, nitrogen deficiency can be noted from chlorosis, the condition in which leaves yellow as the plant's chlorophyll content drops [9]. With reduced photosynthetic activity, the plant will not reach peak health and yield will be low. Water is also key to the transportation of nutrients from the soil to a plant. The availability of water to a plant depends on the weather conditions during the growing season, the soil moisture, and the field micro-topography affecting water flow and accumulation [10,11]. Understanding a field's characteristics as well as monitoring plant biophysical characteristics including height, leaf area, and leaf colour can provide useful information in nitrogen fertilizer applications.
In PA, remote sensing imagery is useful because it does not require physical or destructive contact with plants to gather valuable crop information [12]. Vegetation indices (VIs) can be derived from the spectral information provided by the imager; VIs are mathematical combinations or transformations of spectral bands that have been widely used in agricultural research. VIs allow for the deriving of specific plant properties such as chlorophyll or nutrient content by taking advantage of the differential spectral properties of plants in the visible and near-infrared (NIR) wavelengths [13][14][15]. The VI information can then provide timely knowledge of crop conditions, allowing for a suitable rate of application at the right time and location depending on the variations within a field.
Optical satellite imagery for crop monitoring has had several decades of research and application [16]. Examples of recently launched optical satellites in operation include RapidEye since 2008, Landsat 8 since 2013, and Sentinel-2 since 2015, all frequently used in studies on crop nutrient, yield, and growth management [3]. RapidEye has five spectral bands with 6.5 m resolution. Depending on the location, the five-satellite constellation revisit time is between one to five days. Landsat 8 Operational Land Imager has nine spectral bands with varying spatial resolutions of 15 to 30 m. It has a 16-day revisit time to the same area and takes over 700 scenes a day. Sentinel-2 has 13 spectral bands with 10 m, 20 m, and 60 m spatial resolutions depending on the band. Sentinel-2 constellation is composed of two satellites allowing for a five-day revisit time over the same area. Limitations in optical satellite imagery include low spatial sensitivity as the spatial resolution may be too coarse for small-scale crop fields [12]. The temporal sensitivity can be rather low, such as that of Landsat 8 with a 16-day revisit time, during which crops would have changed significantly, and thus, valuable information on the different stages of growth would not be obtained. Sentinel-2 and RapidEye have higher temporal resolutions of one to five days, but it can vary by location and not all images may be useful due to cloud cover obscuring land. Additionally, factors such as cloud cover, geometric distortion, and atmospheric distortion may require advanced processing expertise to ensure sufficient image quality [17].
New satellite systems are improving in spatial and temporal sensitivity, such as the PlanetScope satellite constellation [18]. Designed for collecting information for use in land-change detection, crop monitoring, climate monitoring, and more, the PlanetScope satellite constellation is composed of over 130 satellites called Doves allowing for spatial resolutions of 3 to 5 m and daily revisit. Beginning with the first launch of a group of Doves in March 2016, over 10 more groups have launched since to improve revisit time, as well as spatial and spectral resolutions. PlanetScope imagery products are also available in multiple asset forms with different radiometric processing and rectification, such as the "surface reflectance" product imagery downloaded for this study. Currently, a select portion of Planet data is available for free download under an open data access policy. PlanetScope imagery has been used in studies of wheat yield, biomass, and LAI monitoring and modelling with promising results [19][20][21]. However, there are few studies focused specifically on nitrogen management using PlanetScope data, a gap that this study aims to fill.
With the rapid advancement in UAV technology in recent years, there is much research interest in UAV-based crop canopy nitrogen retrieval [3]. UAV-based remote sensing can provide low cost and higher spatial and temporal resolution data for crop management. Individuals with basic training can operate a UAV using programmed routes and collect images with <10 cm resolutions [22]. They can be flown to capture more frequent image data and offer flexibility in operation for times when weather is most suitable [23]. Compared to satellites, overall, UAV-based systems are often lower in cost for data collection and Nitrogen 2022, 3 3 processing. Studies have shown significant correlations between crop spectral variables derived from UAV imagery and crop nitrogen content [24][25][26]. Many studies are based on single or combinations of different spectral indices' relationships with crop nitrogen content, noting variation in the relationships at different stages of crop growth [24,27,28]. The spectral indices with the strongest relationship to crop nitrogen were noted to occur during early wheat growth stages before and up to flowering. Often, studies involving the estimation of nitrogen were conducted in controlled experimental conditions, and more studies are needed on real field conditions.
Wheat was selected for this study because it is among the most grown crops in Ontario [16]. With the development of new remote sensing technologies, processing methods, and computing capabilities, estimation models for crop nitrogen can be improved. Machine learning is an area of research interest as it can be used to develop accurate crop monitoring models for large, nonparametric, nonlinear datasets [29]. Recent studies have tested the use of linear regression, Random Forest (RF), and Support Vector Regression (SVR) models in UAV-based canopy nitrogen weight (CNW) prediction models [25][26][27]29]. Although linear regression is a commonly used method to predict nitrogen, some VIs (e.g., NDVI) may saturate beyond the early growth crop stages and some models may have reduced accuracy due to multicollinearity [26,30]. By contrast, machine learningbased regression methods such as RF and SVR were found to produce more accurate models compared to classical linear regression methods, as they are unaffected by the assumptions of linear regression [30]. However, most current literature on the use of remote sensing data and machine learning have only considered spectral information for crop nitrogen modelling [28]. As a crop's nitrogen status can be affected by many factors including fertilizer application, soil characteristics, water availability, and field microtopography, nitrogen prediction models may be improved if these plant physiological and environmental variables are considered [31].
With better management of nitrogen fertilizers, not only can costs and negative environmental impacts be minimized-yield and quality can also be increased. This study aims to evaluate machine learning modelling methods with plant spectral, biophysical, and field environmental variables to predict CNW in wheat crops using UAV and satellite-based imagery. The objectives of this study include, (i) studying the relationship between the spatial variation of CNW and factors such as plant height, LAI, soil moisture, and topographic metrics within wheat fields in Southwestern Ontario using multispectral UAV-or PlanetScope-imagery; (ii) determining the optimal combination(s) of spectral variable(s), crop variables, and/or environmental conditions (soil, water, and topographic data) for wheat canopy nitrogen estimation and prediction; and (iii) evaluating the temporal variation of nitrogen estimation and prediction during the early growth stages of wheat using UAV images or PlanetScope images, and related variables.

Study Area and Data Collection
The study sites are in Strathroy-Caradoc, Ontario, Canada, nearby the community of Mount Brydges (Figure 1). Fieldwork was conducted during May-June 2020 with an average temperature of 22 • C and humidity averaging 73%, characteristic of southwestern Ontario's humid continental climate zone. The sites are in a predominantly agricultural area about 25 km west of London, Ontario's urban center.
In southwestern Ontario, winter wheat is planted in autumn for germination, lying dormant during winter and resuming growth in spring. This study's wheat fields' cultivar was "Soft Red Winter Wheat", which was planted in mid-October 2019, began sprouting mid-April 2020, and was harvested in early August. The three fields labelled W1, W2, and W3 are sized at 48, 21, and 27 hectares, respectively. Beginning in early May, sampling was conducted every 7-8 days for a total of five sample dates to capture significant growth stages of the crop (Table 1). Wheat growth stages were recorded following the Biologische Bundesanstalt, Bundessortenamt, and CHemische Industrie (BBCH) scale [32]. The data Nitrogen 2022, 3 4 acquired for this study encompassed the wheat crop's growth period from leaf development up to and including inflorescence emergence (BBCH 10 to 59) for nitrogen estimation. The early stages of growth before flowering are especially important as fertilizers applied then have better nitrogen use efficiency and yield response [4,33]. Fertilizer can be applied during autumn planting, but lower probability of rainfall also decreases the amount of N from moving into the soil. Fertilizer left above soil during winter will not penetrate and may move during spring snowmelt. During later growth stages (BBCH 60+) from fruit development to ripening, root N uptake slows down as the plants translocate N from vegetation to grains and excess fertilizer can leach from the field [34].
Nitrogen 2021, 2, FOR PEER REVIEW 4 Figure 1. Location of the three 2020 wheat fields near Mount Brydges in southwestern Ontario, Canada, over a Google Earth image.
In southwestern Ontario, winter wheat is planted in autumn for germination, lying dormant during winter and resuming growth in spring. This study's wheat fields' cultivar was "Soft Red Winter Wheat", which was planted in mid-October 2019, began sprouting mid-April 2020, and was harvested in early August. The three fields labelled W1, W2, and W3 are sized at 48, 21, and 27 hectares, respectively. Beginning in early May, sampling was conducted every 7-8 days for a total of five sample dates to capture significant growth stages of the crop (Table 1). Wheat growth stages were recorded following the Biologische Bundesanstalt, Bundessortenamt, and CHemische Industrie (BBCH) scale [32]. The data acquired for this study encompassed the wheat crop's growth period from leaf development up to and including inflorescence emergence (BBCH 10 to 59) for nitrogen estimation. The early stages of growth before flowering are especially important as fertilizers applied then have better nitrogen use efficiency and yield response [4,33]. Fertilizer can be applied during autumn planting, but lower probability of rainfall also decreases the amount of N from moving into the soil. Fertilizer left above soil during winter will not penetrate and may move during spring snowmelt. During later growth stages (BBCH 60+) from fruit development to ripening, root N uptake slows down as the plants translocate N from vegetation to grains and excess fertilizer can leach from the field [34].   Prior to sample point selection, a DJI Phantom 4 Real-Time Kinematics (RTK) UAV was flown over the bare soil of the fields. The UAV connected to an RTK global navigation satellite system (GNSS) base station acquired images with positioning metadata at 1.5 cm vertical and 1 cm horizontal positioning accuracy. These images were mosaicked to create a digital elevation model (DEM) GeoTiff image. The DEM was imported into Google Earth and 16 sample points were selected for each field based on the variation of elevation and coverage of representative areas ( Figure 2). Considering the crop rows were planted northwest-southwest for field W1 and north-south for fields W2 and W3, the sampling pattern followed the row directions for navigation efficiency. The sample points were placed in a four-by-four grid, spaced 60 m from one another. A minimum distance of 50 m from roads was used to reduce possible effects of transportation pollution. From Google Earth, the points were exported to a KML file and downloaded onto mobile devices. On the first fieldwork date, the KML file was used to navigate to the pre-determined sample points and red flags were placed for accurate positioning in following fieldwork sessions.
satellite system (GNSS) base station acquired images with positioning metadata at 1.5 cm vertical and 1 cm horizontal positioning accuracy. These images were mosaicked to create a digital elevation model (DEM) GeoTiff image. The DEM was imported into Google Earth and 16 sample points were selected for each field based on the variation of elevation and coverage of representative areas ( Figure 2). Considering the crop rows were planted northwest-southwest for field W1 and north-south for fields W2 and W3, the sampling pattern followed the row directions for navigation efficiency. The sample points were placed in a four-by-four grid, spaced 60 m from one another. A minimum distance of 50 m from roads was used to reduce possible effects of transportation pollution. From Google Earth, the points were exported to a KML file and downloaded onto mobile devices. On the first fieldwork date, the KML file was used to navigate to the pre-determined sample points and red flags were placed for accurate positioning in following fieldwork sessions. To measure the biomass at each sampling point, a square guide made of plastic tubing of 50 cm × 50 cm was placed around a patch of wheat. The plants within the guide were destructively collected by cutting at the stem base above ground, then placed in paper bags. The average distance between rows was 17 cm, and mostly three rows of wheat would be collected from the 0.25 m 2 biomass block. On the same day following fieldwork collection, the fresh biomass was weighed in grams then fully dried in an oven at 60 °C To measure the biomass at each sampling point, a square guide made of plastic tubing of 50 cm × 50 cm was placed around a patch of wheat. The plants within the guide were destructively collected by cutting at the stem base above ground, then placed in paper bags. The average distance between rows was 17 cm, and mostly three rows of wheat would be collected from the 0.25 m 2 biomass block. On the same day following fieldwork collection, the fresh biomass was weighed in grams then fully dried in an oven at 60 • C for 48 h. Dry biomass was weighed before being sent to A&L Canada Laboratories for plant analysis. The biomass was ground into a fine powder able to pass through a 1 mm sieve before being used in the Laboratory Equipment Company (LECO) FP628 nitrogen combustion method to obtain the leaf N content percentages [35].
Six measurements of plant height in centimeters were taken around each sampling point within a 1 m 2 block on every fieldwork date, and an average height was calculated. Six measurements of soil moisture were collected within a 1 m 2 block around the sampling point using an ML3 ThetaProbe (Delta-T Devices Ltd., Burwell, Cambridge, UK) and averaged [36]. Plant phenology was recorded at each sample point using the BBCH scale to determine growth stage during data collections.
The leaf area index (LAI) was measured non-destructively using a LAI-2200C Plant Canopy Analyzer (Li-Cor, Inc., Lincoln, NE, USA) [37]. Following manufacturer recommen-dations, five measurements were taken along the row at each sample point during clear skies or uniform overcast conditions. Files from the device were transferred to and processed with the Li-Cor software File Viewer FV2200. Processing options include scattering correction based on the field conditions, and output text files with final LAI measurements.

UAV Imagery
In this study, the two UAVs used were a Da-Jiang Innovations (DJI) (DJI, Shenzhen, China) Matrice 100 and a DJI Phantom 4 RTK [38]. The DJI Phantom 4 was flown at 30 m altitude as per the manufacturer's recommendations regarding the visual sensor system's optimal performance and connection to the RTK GNSS base station. The RGB image resolution was 0.9 × 0.9 cm and capture was set to 80% side and 80% front image overlap. A fine resolution, typically < 10 cm, and high image overlap ensure a greater chance of successful image mosaicking as crop canopies densify during a growing season [27].
The DJI Matrice 100 carried a MicaSense RedEdge (MicaSense Inc., Seattle, WA, USA) narrowband multispectral camera including the following bands ordered by increasing wavelength: (1) blue, (2) green, (3), red, (4) red-edge, and (5) near-infrared (NIR) ( Table 2) [39]. Image acquisition was conducted on the same dates as ground fieldwork before biomass collection. For the fieldwork conducted during the week of 27 May, the weather conditions were characterized by strong winds and sudden rain showers, so the UAV was flown on different dates between 26 May and 29 May for each field whenever the weather was suitable. The flight plans for each field were made in the Pix4D software suite "Pix4Dcapture" app (Pix4D S.A., Prilly, Switzerland) to cover entire fields in a zigzag pattern [40]. Pix4Dcapture includes the function of adding custom UAV and sensor properties to calculate the flight plan's estimated total time, battery usage, and image resolution at the plan's designated altitude. The flight altitude for each field's plan was set at 60 m, as >60 m wind and gust conditions were usually greater than the UAV manufacturer recommendations. Based on the MicaSense camera's specifications, the resulting image resolution was 4 cm × 4 cm with capture set to 80% side and 80% front overlap. Unfortunately, images for field W3 on 5 May and field W1 on 27 May were corrupted and unusable for further processing.

UAV Image Processing
UAV image processing followed the flowchart of Figure 3. The DJI Phantom 4 RTK images taken in April over each field's bare soil were mosaicked in Pix4Dmapper software to generate a DEM GeoTiff. Using the System for Automated Geoscientific Analysis (SAGA), a free, open-source spatial data analysis software package, topographic metrics were computed from the DEM [41]. Using the "Terrain Analysis-Morphometry" tool, the metrics generated include slope, aspect, profile curvature, and plan curvature. Creating topographic wetness indices (TWI) required several steps. First, using the "Terrain Analysis-Hydrology" tool, a flow accumulation layer of the field was created. Using the "Flow Width and Specific Catchment Areas (SCA)" tool, two SCAs were created using different algorithms: Deterministic 8 and Multiple Flow Direction. The "Topographic Wetness Index (TWI)" tool was then used to create two maps, TWI #1 and TWI #2, based on the respective algorithm SCAs. The final topographic metrics were exported as GeoTiff files.
ing topographic wetness indices (TWI) required several steps. First, using the "Terrain Analysis-Hydrology" tool, a flow accumulation layer of the field was created. Using the "Flow Width and Specific Catchment Areas (SCA)" tool, two SCAs were created using different algorithms: Deterministic 8 and Multiple Flow Direction. The "Topographic Wetness Index (TWI)" tool was then used to create two maps, TWI #1 and TWI #2, based on the respective algorithm SCAs. The final topographic metrics were exported as GeoTiff files. MicaSense camera multispectral images were processed in Pix4Dmapper to create an orthomosaic image per band with 4 cm × 4 cm resolution. An important step in producing a high-quality final image is radiometric calibration considering the sensor influence and scene illumination of the UAV flight. Prior to each flight over a field, the MicaSense camera was positioned above a MicaSense Calibrated Reflectance Panel to take a minimum of five white reference images for each band. From the Exchangeable Image File Format (EXIF) metadata of the images, Pix4Dmapper can read the sensor settings, properties, and geolocation at the time the images were taken. Before starting the mosaic processing, setting options include inputting the white reference images and manufacturer-provided panel reflectance values to calibrate and correct image reflectance for each of the bands. Then, Pix4Dmapper uses the Structure from Motion (SfM) algorithm to correct image perspectives and recognize where to stitch images together [22]. The high image overlap parameters set in the flight plan enable the software to recognize greater, similar areas of MicaSense camera multispectral images were processed in Pix4Dmapper to create an orthomosaic image per band with 4 cm × 4 cm resolution. An important step in producing a high-quality final image is radiometric calibration considering the sensor influence and scene illumination of the UAV flight. Prior to each flight over a field, the MicaSense camera was positioned above a MicaSense Calibrated Reflectance Panel to take a minimum of five white reference images for each band. From the Exchangeable Image File Format (EXIF) metadata of the images, Pix4Dmapper can read the sensor settings, properties, and geolocation at the time the images were taken. Before starting the mosaic processing, setting options include inputting the white reference images and manufacturer-provided panel reflectance values to calibrate and correct image reflectance for each of the bands. Then, Pix4Dmapper uses the Structure from Motion (SfM) algorithm to correct image perspectives and recognize where to stitch images together [22]. The high image overlap parameters set in the flight plan enable the software to recognize greater, similar areas of each image for a higher chance of mosaic success. Pix4Dmapper processing options include georeferencing with ground control points (GCP) to improve the absolute location accuracy. Following Pix4D's recommended number of five GCPs, black and white checkered boards were placed around the wheat fields and their coordinates recorded using a Global Positioning System (GPS) connected to the RTK. The output orthomosaic images are GeoTiff files with reflectance values for each MicaSense band.

Satellite Imagery
In this study, PlanetScope satellite imagery was acquired from Planet Labs Inc. through submission and acceptance of the project proposal. Revisit times are almost daily worldwide with resampled spatial resolution of 3 m × 3 m. The image products available for this study were from the third-generation sensors, Dove-R, with four bands ordered by increasing wavelength: (1) blue, (2) green, (3) red, (4) NIR (Table 3). The product scenes were approximately 25 km × 23 km , and in the case of the three study fields, all were in one image. Five satellite images were downloaded, matching the capture dates with the ground data collection dates to maintain consistency with the plant physiology and field conditions. The images were available for download as GeoTiff surface reflectance assets, orthorectified and radiometrically corrected based on the atmospheric conditions of the Nitrogen 2022, 3 8 specific ground locations. In ArcGIS, the large scenes were cropped to smaller images of each study field for subsequent processing.

Vegetation Indices
MicaSense and PlanetScope reflectance images of each of their respective bands were exported into ArcGIS to extract crop canopy reflectance values at the sample points. With the "Raster Calculator" tool in PCI Geomatica Banff, the images were used to compute the VIs listed in Table 4 and exported as GeoTiff files. Using ArcGIS, the VI layer values were extracted at each sample point. As MicaSense has a red-edge band, the following three indices were included specifically for the UAV-data modelling: the chlorophyll index rededge (CI_RE), the normalized difference vegetation index (NDRE), and the ratio vegetation index #2 (RVI2). VIs have been extensively studied for the purpose of crop monitoring and biophysical estimation, such as the normalized difference vegetation index (NDVI) [12,42]. VIs developed for chlorophyll estimation have been found to be related to plant nitrogen content as the photosynthetic enzyme, rubisco, encompasses the largest proportion of nitrogen in leaves [4]. Chlorophyll reflects green and NIR radiation and absorbs more than 70% of blue and red radiation [43].

Canopy Nitrogen Weight Modelling
CNW was calculated using the following method [58]: Nitrogen 2022, 3 9 where CNW is the canopy nitrogen weight (g/m 2 ), LNC is the leaf nitrogen content (%), and W d is the dry biomass weight (g/m 2 ). CNW assumes the plants collected for biomass around the sample point, within the 0.25 m 2 block, have the same amount of nitrogen. The total biomass per area was used as the dry biomass of wheat plants at early growth stages (BBCH < 60) since the leaves constitute most of the plant weight. Compared to other biophysical parameters such as plant nitrate content and plant nitrogen concentration (%), CNW has been found to have greater correlation with spectral data [28]. RF is an ensemble learning method suited for the classification or regression of large, nonparametric datasets. Training data are randomly selected from the dataset based on a user-defined percentage; commonly, 70% of a dataset is used to train a model. From the training data, many decision trees are generated by the algorithm to determine the variables' importance in the regression. Decision trees split at nodes depending on the independent variable that contributes most to the dependent variable. The validation data are the remainder of the dataset that is not used in training, and the average output from decision trees is used to evaluate the model's performance. Advantages of RF include quick computation, no overfitting from training data, and high performance in studies [29]. Support Vector Machines (SVM) are supervised learning algorithms used for classification and regression. SVR uses a decision boundary, known as a hyperplane, to split classes of training data based on data characteristics. The data points closest to either side of the hyperplane are known as support vectors, used as training samples to determine the optimal hyperplane position from the midpoint of the margin. SVR performs modelling in a high-dimensional space. SVR uses a kernel trick (i.e., the Radial Basis Kernel) for nonlinear data, placing the data in a dimensional space to separate into groups based on radial distance between data points. SVR has advantages over simple linear regression models, as its flexibility with nonparametric data has better modelling capabilities [59]. The modelling for this study was written in the R programming language using R Studio (Version 4.0.3), an open-source and free Integrated Development Environment (IDE) [60]. For RF, the "randomForest" package was used, and for SVR, the "e1071" package was used.
Modelling was performed separately for MicaSense-based and PlanetScope-based sensor variables for comparison of UAV and PlanetScope data results. The MicaSense models used all five MicaSense bands, all 14 VIs listed in Table 1, plant height, LAI, soil moisture, and topographic metrics, and the dependent variable was the CNW. The PlanetScope models used all four PlanetScope bands, 11 VIs listed in Table 1 not including the red-edge based indices (CI_RE, NDRE, RVI2), plant height, LAI, soil moisture, and topographic metrics, and the dependent variable was the CNW. Measurements of variables at each sample point were averaged to a 1 m 2 scale. Datasets were randomly divided into 70% training set and 30% validation set. The quality of the models was assessed using the coefficient of determination (R 2 ) and the Root Mean Square Error (RMSE), calculated using Equations (2) and (3), respectively [61]: where y i is the observed value,ŷ i is the predicted value, and y is the mean of the observed values of the dataset; and: whereŷ i is the predicted CNW value (g/m 2 ), y i is the observed CNW value (g/m 2 ), n is the number of observations, and i is the index of summation in increments of one.

Nitrogen Statistics
Overall, the canopy N weight for the wheat fields increased in variation during the fieldwork season (Figure 4). There was a slight decrease in canopy N weight from 6 to 12 May, likely due to a period of several days of continuous rainfall between samplings. Rainfall can lead to leaching of soil N from areas around plant roots, in this case reducing plant N utilization while the biomass continued to increase [62,63]. After 12 May, the farmer applied fertilizer to the wheat field once before 20 May and, consequently, the canopy N weight increased during the following weeks. There was an outlier from a sample point on 20 May, which was removed from the dataset for modelling.
Nitrogen 2021, 2, FOR PEER REVIEW 11 Figure 4. Variation of canopy nitrogen weight (g/m 2 ) as a function of the date of field measurements during the 2020 growing season. An outlier is represented by a dot on the graph.

UAV Regression Models
For the UAV regression models, the first modelling step used all 28 variables including: the five MicaSense band reflectances, 14 VIs, plant height, soil moisture, LAI, topographic slope, aspect, profile curvature, plan curvature, TWI #1, and TWI #2. Single-date and combinations of multi-date datasets were tested to evaluate the temporal effect on models. A proportion of 70% of each dataset was used to calibrate the RF and SVR models, before the remaining 30% was used for validation. From the calibrated models (

Soil Statistics
A VitTellus soil health analysis [64] was conducted on the study fields, producing the following mean values: soil nitrate-nitrogen of 46.81 ppm, mineralizable nitrogen of 29.17 ppm, and water-extracted soil nitrate of 44.38 ppm. A&L Canada Soil Health Index ratings were in the "Good-High" category, and the soil textural class for the fields was predominantly sandy loam.

UAV Regression Models
For the UAV regression models, the first modelling step used all 28 variables including: the five MicaSense band reflectances, 14 VIs, plant height, soil moisture, LAI, topographic slope, aspect, profile curvature, plan curvature, TWI #1, and TWI #2. Single-date and combinations of multi-date datasets were tested to evaluate the temporal effect on models. A proportion of 70% of each dataset was used to calibrate the RF and SVR models, before the remaining 30% was used for validation. From the calibrated models (Table 5), RF had better performance in comparison to SVR; specifically, it had higher R 2 and lower RMSE values overall. The top model performance was that of RF with the combination of 12, 20, and 27 May resulting in an R 2 of 0.96 and an RMSE of 1.07 g/m 2 . Although the 20/27 May and 4 June model had a slightly higher R 2 of 0.97, the RMSE was considerably greater at 1.76 g/m 2 . Of the RF and SVR models, that of 20 May had the lowest performance. Of the SVR models, 12, 20, and 27 May model had the best performance with an R 2 of 0.87 and an RMSE of 2.07 g/m 2 . The UAV models were then applied to the validation datasets ( Table 6). All single date models for RF and SVR had low performance. This was likely due to the small number of sample point data from single date sets resulting in the calibration model not encompassing normal field data variation. The top performing model was RF with the 12, 20, and 27 May data, resulting in an R 2 of 0.74 and an RMSE of 2.76 g/m 2 . Overall, RF had better performance than SVR and multi-date models had better performance than single-date models.

PlanetScope Regression Models
For the PlanetScope regression models, 24 variables were used in the first step, including: the four PlanetScope band reflectances, 11 VIs, plant height, soil moisture, LAI, topographic slop, aspect, profile curvature, plan curvature, TWI #1, and TWI #2. To test the temporal effects on the models, single-date and multi-date dataset combinations were used. Taking 70% of the dataset to calibrate the RF and SVR models, overall RF had better performance to SVR with higher R 2 and lower RMSE values (Table 7). Of all models, RF with the combination of the 12, 20, and 27 May data had best performance with an R 2 of 0.96 and an RMSE of 1.10 g/m 2 . Of the SVR models alone, the 12, 20, and 27 May model had the highest performance with an R 2 of 0.87 and an RMSE of 2.12 g/m 2 . Next, the PlanetScope models were applied to the remaining 30% of datasets for validation (Table 8). Of all the models, single-date 20 May and 4 June had very low performance for both RF and SVR. Comparing all RF and SVR models, RF usually had better performance than SVR except in the single date models of 5 May, 12 May, and 20 May, and the multi-date model with the three dates combined. The top performing model was RF with 12, 20, and 27 May data resulting in an R 2 of 0.83 and an RMSE of 1.77 g/m 2 .

Variable Importance Plots
The top performing UAV and PlanetScope regression models with all variables were both from RF and were multi-date combinations of the 12, 20, and 27 May datasets. Modelling of RF in RStudio can be visualized with a variable importance plot using the "varImp-Plot()" function. The higher a variable's IncNodePurity value is, the more important the explanatory variable is in the prediction of CNW. In the UAV model plot ( Figure 5) with 28 variables, plant height was by far the most important predictor, followed by LAI. Rededge VIs (NDRE, RVI2, CI_RE) were the third, fourth, and sixth most important predictor variables, respectively. Of the MicaSense band reflectances, NIR was the most important. Soil moisture also appeared to be among the top group of important variables with a greater difference in IncNodePurity compared to the other variables in the plot below. Of the topographic metrics, TWI #2 and TWI #1 held more importance than the rest.  Table 1 for the full name of vegetation indices. TWI_1, total wetness index #1; TWI_2, total wetness index #2.
The satellite RF regression model plot with 24 variables indicates height as the most important predictor for CNW ( Figure 6). Of the four PlanetScope band reflectances, the blue band was most important and second-most important overall of all variables, followed closely by LAI. Overall, VIs and PlanetScope bands were among the most important variables, while topographic metrics and soil moisture were the least important.  Table 1 for the full name of vegetation indices. TWI_1, total wetness index #1; TWI_2, total wetness index #2.
The satellite RF regression model plot with 24 variables indicates height as the most important predictor for CNW ( Figure 6). Of the four PlanetScope band reflectances, the blue band was most important and second-most important overall of all variables, followed closely by LAI. Overall, VIs and PlanetScope bands were among the most important variables, while topographic metrics and soil moisture were the least important.  Table 1 for the full name of vegetation indices. TWI_1, total wetness index #1; TWI_2, total wetness index #2.

UAV Regression Models
From the variable importance plot of the best-performing UAV regression model using 12, 20, and 27 May data, combinations of variables were tested in the RF and SVR models. Evaluating based on thresholds of variable importance, as shown in Figure 5, we selected variables for the top 7, 8, 9, 13, and 16 variable groups. An additional group containing only MicaSense band reflectances and VIs was included for comparison as common modelling approaches in other studies use only spectral variables [3]. Furthermore, 70% of each variable group dataset was used for model calibration (Table 9). Overall, RF had better performance than SVR in terms of higher R 2 and lower RMSE values. All variable group models using RF had a high R 2 of 0.96 except for the spectral-only variable group with lower R 2 . RF models' RMSE values were also very close, within <0.10 g/m 2 difference, the lowest being the top 9 and top 13 variable groups at 1.07 g/m 2 . Again, the RF spectral-only variable group had poorer performance in comparison to other RF models, as its RMSE value was significantly higher. The SVR spectral-only variable group performed poorly compared to all the other SVR models. Of the SVR models, the top 7 variable group had the best performance with an R 2 of 0.86 and an RMSE of 2.05 g/m 2 . The top 8,9,13, and 16 variable SVR models had only slightly lower R 2 and higher RMSE values compared to the best SVR model.  Table 1 for the full name of vegetation indices. TWI_1, total wetness index #1; TWI_2, total wetness index #2.

UAV Regression Models
From the variable importance plot of the best-performing UAV regression model using 12, 20, and 27 May data, combinations of variables were tested in the RF and SVR models. Evaluating based on thresholds of variable importance, as shown in Figure 5, we selected variables for the top 7, 8, 9, 13, and 16 variable groups. An additional group containing only MicaSense band reflectances and VIs was included for comparison as common modelling approaches in other studies use only spectral variables [3]. Furthermore, 70% of each variable group dataset was used for model calibration (Table 9). Overall, RF had better performance than SVR in terms of higher R 2 and lower RMSE values. All variable group models using RF had a high R 2 of 0.96 except for the spectral-only variable group with lower R 2 . RF models' RMSE values were also very close, within <0.10 g/m 2 difference, the lowest being the top 9 and top 13 variable groups at 1.07 g/m 2 . Again, the RF spectral-only variable group had poorer performance in comparison to other RF models, as its RMSE value was significantly higher. The SVR spectral-only variable group performed poorly compared to all the other SVR models. Of the SVR models, the top 7 variable group had the best performance with an R 2 of 0.86 and an RMSE of 2.05 g/m 2 . The top 8, 9, 13, and 16 variable SVR models had only slightly lower R 2 and higher RMSE values compared to the best SVR model.
The calibrated models were applied to the remaining 30% of datasets for validation (Table 10). In the validation models, SVR had better performance than RF except for the spectral-only model. The best-performing model was SVR with the top seven variables, resulting in an R 2 of 0.80 and an RMSE of 2.62 g/m 2 . Compared to the best model, other SVR models had close performance, but as more variables were added, the R 2 had lower values.

PlanetScope Regression Models
For the UAV regression models using single and multi-date datasets, the best-performing combination was the model using the 12, 20, and 27 May data with RF. Evaluating based on thresholds of variable importance, as shown in Figure 6, the selected variable groups for further model testing included the top 6, 10, 13, 17, and a group containing only PlanetScope band reflectances and VIs. A proportion of 70% of each variable group dataset was used for model calibration (Table 11). Of the calibration models, all RF models had the same high R 2 value with low RMSE compared to all SVR models. The best model performance was RF with the top six variables, with an R 2 at 0.96 and an RMSE 1.18 g/m 2 . For the SVR models, the best-performing was the combination of the top six variables, resulting in an R 2 of 0.84 and an RMSE of 2.31 g/m 2 . From the calibration models, the spectral-only SVR model had the poorest performance with an R 2 of 0.77 and an RMSE of 2.81 g/m 2 . The remaining 30% of variable datasets were used for the validation models (Table 12). The poorest model performances for RF and SVR were from the spectral-only variable group, resulting in the lowest R 2 and highest RMSE values out of all tested. SVR had better performance compared to RF in the spectral-only, top 6, and top 10 groups, but for the top 13 and top 17 groups, RF was better. The best-performing model was RF with the top 17 variables, with an R 2 of 0.92 and an RMSE of 1.75 g/m 2 .

Crop Nitrogen Prediction Maps
To create the crop nitrogen prediction maps, raster layers of each variable were required. The VI raster layers were created using the raster calculator in PCI Geomatica Banff from the MicaSense band rasters mosaicked from Pix4Dmapper. To extract raster layers for height, the Phantom 4 RTK flight processed in Pix4Dmapper includes the option for a digital surface model (DSM) output. The DSM captures the natural and built features of the environment. By subtracting the DEM raster from the DSM raster pixel values using ArcGIS raster calculator, the output raster layer has the height data of the crop [65].
Generating an LAI raster layer is a more involved process. The process is based on methodology proposed by Song et al. (2020) using a simulated observation of a point cloud Multiview angle (SOPC-M) designed to obtain a 3D spatial distribution of vegetation and bare ground points to calculate the gap fraction, followed by the crop's effective LAI (LAIe), from a UAV-based 3D point cloud dataset [66]. In ArcGIS, the point data layer was converted to a raster layer using the "Point to Raster" conversion tool to a 2 m × 2 m resolution. From the linear regression equation, the digital hemispherical photography (DHP) LAIe could be calculated using the ArcGIS raster layer to create the final LAI raster used in the model. For field W3, the LAI layer processing area error unfortunately only covered 12 of the 16 sample points; thus, the final map produced did not cover all sample points.
Prior to running the regression models in R, the raster layers needed to have the same resolution and extent for the functions to stack them. All raster layers were run through the ArcGIS "Resample" tool to 1 m × 1 m resolution, as the CNW variable is based on a 1 m 2 area. Resampling was conducted with the bilinear technique, which is a bilinear interpolation and determines the new value of a cell based on a weighted distance average of the four nearest input cell centers. As the LAI raster extent was the smallest of all rasters, it was used as the output extent reference layer for the ArcGIS "Clip Raster" tool with the selected option to maintain clipping extent. The columns and rows of the output raster were adjusted, and the pixels were resampled to match the reference layer exactly. Next, the ArcGIS "Extract by Mask" tool was used to extract the cell values specifically where the LAI raster extent pixels held values greater than zero. Thus, each variable layer was prepared for the R modelling.
For the UAV best-performing model, SVR with the top seven variables included height, LAI, NDRE, RVI2, BNDVI, CI_RE, and NIR. In R, the "raster", "rgdal", and "rasterVis" libraries were used to generate the final prediction map outputs. The variable raster layers were stacked using the "raster::stack" function, and the "raster::predict" tool was used with the selected UAV-based SVR model using the top seven variables. The resulting prediction rasters for each field were exported as GeoTiff files into ArcGIS to create the final prediction maps with 1 m × 1 m resolution (Figure 7). The low and high CNW areas are displayed in red and green, respectively, for distinct contrast between the nitrogen levels. From the prediction rasters, the predicted CNW values around each sample point were extracted and compared to the measured values in-field. The resulting RMSE values for W1, W2, and W3 were 4.27, 2.32, and 3.08 g/m 2 , respectively.
The best-performing satellite model used RF with the top 17 variables: height, LAI, 11 VIs, and the four PlanetScope bands. Following the same processing steps as the UAVbased model, the satellite variable raster images were processed to have the same extent and 1 m × 1 m resolution. As the satellite images began with 3 m × 3 m resolution, much larger compared to the UAV images, the resulting raster layers had a smoothed appearance. The final prediction rasters for each field are displayed as maps in Figure 8 with 1 m × 1 m resolution. From the prediction rasters, the predicted CNW values around each sample point were extracted and compared to the measured values in-field. The resulting RMSE values for W1, W2, and W3 were 3.12, 1.79, and 3.08 g/m 2 , respectively.
Comparing the maps in Figures 7 and 8, it can be seen that there were differences in the spatial distribution and values of CNW. It appears that the UAV maps show more areas with low canopy nitrogen weight, while PlanetScope shows more areas with high CNW. This is likely due to the resampling method, as the UAV's centimeter-level image resolution includes spectral data of the visible soil between crop rows, while for PlanetScope images with coarser meter-level resolution, the pixel values average towards the spectral values of the crop rather than the soil. The UAV-based nitrogen prediction maps may have lower performance compared to satellite-based maps because of the volatile weather conditions during the fieldwork conducted on the week of May 27th. The UAV was flown on different dates between May 26th and May 29th for each field whenever the weather was suitable, and the crop spectral characteristics may have changed. However, the satellite images for the three fields for May 27th were all from the same date and present a more consistent relationship with the fieldwork data. Nitrogen 2021, 2, FOR PEER REVIEW 20 Comparing the maps in Figures 7 and 8, it can be seen that there were differences in the spatial distribution and values of CNW. It appears that the UAV maps show more areas with low canopy nitrogen weight, while PlanetScope shows more areas with high CNW. This is likely due to the resampling method, as the UAV's centimeter-level image resolution includes spectral data of the visible soil between crop rows, while for Plan-etScope images with coarser meter-level resolution, the pixel values average towards the spectral values of the crop rather than the soil. The UAV-based nitrogen prediction maps may have lower performance compared to satellite-based maps because of the volatile weather conditions during the fieldwork conducted on the week of May 27th. The UAV was flown on different dates between May 26th and May 29th for each field whenever the weather was suitable, and the crop spectral characteristics may have changed. However, the satellite images for the three fields for May 27th were all from the same date and present a more consistent relationship with the fieldwork data.

Discussion
For this study, the RF and SVM regression methods were used to predict the CNW of wheat using UAV MicaSense band reflectances, PlanetScope band reflectances, selected VIs, plant height, LAI, soil moisture, and topographic metrics. The models created were grouped according to UAV-based and satellite-based data.
For the UAV RF and SVR regression models, calibration was conducted with 28 variables from single and multi-date datasets. Evaluating the validation models of each dataset, the performance of UAV single-date models was poor with R 2 values of, at most, 0.25 and overall non-significant results. Combining UAV multi-date data yielded better results, with the best performance from the RF three-date model of 12, 20, and 27 May resulting in an R 2 of 0.74 and an RMSE of 2.76 g/m 2 . For the PlanetScope RF and SVR models, the calibration of the models used 24 variables for single and multi-date datasets. Of the validation models, the single-date models of 20 May and 4 June had the lowest performance. However, the other PlanetScope single-date models' performances were much better overall compared to the UAV single-date models. In general, the PlanetScope multi-date models did not have significantly better results than its single-date models. The best-performing PlanetScope model was that which was based on three dates, 12, 20, and 27 May, using RF, with an R 2 of 0.83 and an RMSE of 1.77 g/m 2 . Both the best UAV and best satellite models were from 12, 20, and 27 May data, during which the wheat crops were in the BBCH 23-41 growth stages, mainly defined by tillering, stem elongation, and

Discussion
For this study, the RF and SVM regression methods were used to predict the CNW of wheat using UAV MicaSense band reflectances, PlanetScope band reflectances, selected VIs, plant height, LAI, soil moisture, and topographic metrics. The models created were grouped according to UAV-based and satellite-based data.
For the UAV RF and SVR regression models, calibration was conducted with 28 variables from single and multi-date datasets. Evaluating the validation models of each dataset, the performance of UAV single-date models was poor with R 2 values of, at most, 0.25 and overall non-significant results. Combining UAV multi-date data yielded better results, with the best performance from the RF three-date model of 12, 20, and 27 May resulting in an R 2 of 0.74 and an RMSE of 2.76 g/m 2 . For the PlanetScope RF and SVR models, the calibration of the models used 24 variables for single and multi-date datasets. Of the validation models, the single-date models of 20 May and 4 June had the lowest performance. However, the other PlanetScope single-date models' performances were much better overall compared to the UAV single-date models. In general, the PlanetScope multi-date models did not have significantly better results than its single-date models. The best-performing PlanetScope model was that which was based on three dates, 12, 20, and 27 May, using RF, with an R 2 of 0.83 and an RMSE of 1.77 g/m 2 . Both the best UAV and best satellite models were from 12, 20, and 27 May data, during which the wheat crops were in the BBCH 23-41 growth stages, mainly defined by tillering, stem elongation, and the beginning of the booting stage. As noted by Hawkesford (2017), the application of nitrogen fertilizer during these early growth stages before flowering is most conducive to efficient nitrogen use and yield response [67]. The ability to accurately estimate the nitrogen levels of crops during early growth stages would be most beneficial for farmers.
In the RF variance importance plot of the best-performing UAV model, of all variables, plant height was the most important predictor of CNW. Song and Wang (2019) also noted that plant height is useful in estimating phenology, biomass, and yield in addition to nitrogen uptake in wheat [68]. On the plot, LAI was the second most important predictor of CNW. LAI has been used extensively in studies to successfully predict crop chlorophyll content, biomass, and yield [69,70]. The study by Zhao et al. (2014) found a significant positive relationship between LAI and differences in crop nitrogen content across wheat growth stages [71]. The gap fraction method of calculating LAI is more accurate during the earlier growth stages of a crop when the canopy is not so dense, allowing for contrast between the vegetation and soil or vegetation and sky images [68]. Among the VI's used in the model, the red-edge VIs (NDRE, RVI2, and CI_RE) were amongst the most important. The red-edge region (680-800 nm) has been shown to encompass sharp changes in the canopy reflectance and can be used to identify important biophysical parameters of the crop. Nitrogen levels have shown the sensitivity of the red-edge region in estimating leaf chlorophyll content due to the high absorption of red radiation and high reflectance of NIR radiation [69,72]. Of the MicaSense bands individually, the NIR band was of highest importance in the model while other individual bands had little effect. Soil moisture also appeared as a variable of high importance, and subsequent variables on the plot had noticeably lower importance. Studies have noted the importance of soil moisture in soil nitrogen mineralization, crop nitrogen uptake, and utilization [73,74]. Of the topographic metrics, the topographic wetness indices were most important, while the remaining metrics had little effect.
In the RF variance importance plot of the best-performing satellite model, similar to the UAV model, plant height was the most important predictor of CNW, with LAI following closely. Interestingly, the PlanetScope blue band was the second most important variable. As the PlanetScope blue band has greater width compared to MicaSense, perhaps the wider bandwidth captured a change of canopy reflectance in the blue-edge region (480-517 nm) which was previously noted in the study by Wei et al. (2008) to be related to crop nitrogen [75]. Other PlanetScope bands in the model had varying levels of importance interspersed amongst the 11 VIs used. On the plot, other non-spectral variables of soil moisture and topographic metrics were of least importance to the model.
From the UAV-based RF variance importance plot, groups of variables were selected for testing in models. Groups of the top 7, 8, 9 13, 16, and spectral-only variables were modeled, with the group of top seven variables of the SVR model performing the best with an R 2 of 0.80 and an RMSE of 2.62 g/m 2 . The top seven variables included plant height, LAI, all three red-edge VIs, BNDVI, and the MicaSense NIR band. Compared to the best models from studies by Asataoui et al.  [24,26,27]. With the UAV-based best model in this study, significantly lower RMSE is a major advantage in terms of reducing the costs of nitrogen fertilizer.
The satellite-based variance importance plot was used to select variable groups for model testing including the top 6, 10, 13, 17, and spectral-only variables. The RF model group of the top 17 variables had the best performance with an R 2 of 0.92 and an RMSE of 1.75 g/m 2 . The height, LAI, all four PlanetScope bands, and total 11 VIs were the variables in the best-performing model. For both the UAV and satellite best-performing selected-variable models, plant height and LAI were the only non-spectral variables. With methods for deriving the height and LAI of a wheat crop field from the UAV Phantom 4 RTK imagery, all variables in the top models can be obtained from in-situ, non-destructive, remote sensing data.
For both the UAV and satellite spectral-only variable groups, the results were poor with R 2 values <0.50 and significantly higher RMSE values compared to other tested variable groups. This is consistent with the studies by Astaoui et al. (2021) and Schirrmann et al. (2016), in which it was noted that, within wheat crops, UAV imagery was limited for the observation of nitrogen status but had good performance in the monitoring of biophysical parameters [27,28]. For the UAV validation spectral-only model, the RMSE dropped by 32% as compared to the best-performing top seven variable model. In the satellite validation spectral-only model, the RMSE dropped by 45% compared to the best-performing top 17 variable model.
In the final validation of canopy nitrogen models with variable combinations, the UAV SVR models mostly had greater R 2 values but also greater RMSE values compared to the RF models. In the UAV spectral-only variable group models, RF had better results than SVR. Considering studies with spectral-only variables for crop nitrogen models, the results are consistent, with RF yielding better nitrogen level prediction compared to SVR models [25,29,76]. Only in the UAV best-performing model of top seven variables was SVR performance better in terms of having both a higher R 2 value and a lower RMSE compared to RF. Of the satellite variable combination models, SVR had better performance than RF except for the top 13 and 17 variable groups. The best-performing satellite model was RF with the top 17 variable group. Although it appears difficult to determine if RF or SVR models are better when built with non-spectral and spectral variables together, ultimately the ideal result is a model which can most accurately predict canopy nitrogen in wheat. In both the UAV and satellite models with different variable combinations, the top variable groups had good overall performances. In comparison to studies with spectral-only variable models, the variable combination models in this study all had lower RMSE values. In the context of nitrogen estimation and practical application, lower RMSE (g/m 2 ) in models is most beneficial for fertilizer management recommendations.

Conclusions
In this study, machine learning regression methods were tested to predict wheat CNW using UAV MicaSense band reflectances, PlanetScope band reflectances, associated VIs, plant height, LAI, soil moisture, and topographic metrics. For UAV models using 28 variables, the combination of 12, 20, and 27 May data with the RF validation model produced the best results with an R 2 of 0.74 and an RMSE of 2.76 g/m 2 . From the model's variable importance plot, the top 7, 8,9,13,16, and spectral-only variable groups were tested. The best validation model used SVR with the top seven variables, which included plant height, LAI, four VIs, and the MicaSense NIR band. For the PlanetScope models using 24 variables, the best performing model was RF with 12, 20, and 27 May data resulting in an R 2 of 0.83 and an RMSE of 1.77 g/m 2 . Based on the model's variable importance plot, the top 6, 10, 13, 17, and spectral-only variable groups were tested. The validation model with the best performance was RF using the top 17 variables including height, LAI, all four PlanetScope bands, and 11 VIs.
A common limitation of in-situ agricultural models, including those developed in this study, are their empirical nature, and applicability can be limited to the dataset they are built and validated upon. Each field and growing season has different conditions and factors that affect plant growth, so models will need further testing to determine their effectiveness in precision agriculture methods. PlanetScope satellite constellation has also launched, and a third generation of sensors were introduced in 2020, known as SuperDove, with the potential to capture imagery with eight spectral bands including a red-edge band. Future work can consider further testing satellite-based nitrogen prediction models including red-edge VI variables.