Estimating Rangeland Forage Production Using Remote Sensing Data from a Small Unmanned Aerial System (sUAS) and PlanetScope Satellite

: Rangelands cover ~23 million hectares and support a $3.4 billion annual cattle industry in California. Large variations in forage production from year to year and across the landscape make grazing management difﬁcult. We here developed optimized methods to map high-resolution forage production using multispectral remote sensing imagery. We conducted monthly ﬂights using a Small Unmanned Aerial System (sUAS) in 2017 and 2018 over a 10-ha deferred grazing rangeland. Daily maps of NDVI at 30-cm resolution were ﬁrst derived by fusing monthly 30-cm sUAS imagery and more frequent 3-m PlanetScope satellite observations. We estimated aboveground net primary production as a product of absorbed photosynthetically active radiation (APAR) derived from NDVI and light use efﬁciency (LUE), optimized as a function of topography and climate stressors. The estimated forage production agreed well with ﬁeld measurements having a R 2 of 0.80 and RMSE of 542 kg/ha. Cumulative NDVI and APAR were less correlated with measured biomass ( R 2 = 0.68). Daily forage production maps captured similar seasonal and spatial patterns compared to ﬁeld-based biomass measurements. Our study demonstrated the utility of aerial and satellite remote sensing technology in supporting adaptive rangeland management, especially during an era of climatic extremes, by providing spatially explicit and near-real-time forage production estimates. study area during 2017-2018 growing season. The time series shows a different phenology for plants growing on south- and north-facing slopes.


Introduction
Rangelands are a key global resource, both in terms of extent and ecological and economic impact. In California, over 60% of the land area is rangeland, which provides a wide range of ecological services, including forage production for livestock and wildlife, water quality protection, biodiversity, recreation, and wildlife habitat [1]. More than half of California rangelands are grazed [2], supporting a $3.4-billion per year livestock industry [3]. Rangeland forage production is characteristically variable and depends on a multitude of drivers-most notably climate, soils and topography, which together regulate plant-available water [4]. Most California rangelands are rain fed, and thus highly vulnerable to irregular precipitation patterns [5]. For example, large year-to-year variations in forage production (from 33 to 11503 kg/ha) were observed in a 14-year record from the California Central Coast region [6]. Warming, drought, and increasing climate variability [7,8] are predicted to become more prevalent this critical metric needs to be calibrated rigorously in different ecosystems and with different input datasets, because of its high model sensitivity.
Remotely sensed imagery from satellites provides a relatively reliable source for obtaining spatiotemporal inputs for LUE-based models. Landsat remote sensing data, with a 30-m resolution and a 16-day revisiting frequency, have been widely used for many vegetation monitoring studies worldwide [23,[33][34][35][36]. However, coincidence of rainfall and growing season in California's Mediterranean climate greatly limits the availability of cloud-free Landsat images during critical plant growth stages from November to April. Moreover, the spatial resolution of traditional satellite remote sensing generally ranges from a kilometer to 10s of meters, whereas most field biomass measurements are at a scale of 10s of centimeters. This scale mismatch poses challenges for developing and validating the LUE optimization.
Recent advances in small unmanned aerial system (sUAS) technology and image processing make it possible to overcome many challenges involved in quantifying forage production across the landscape. Compared to satellites, sUAS can be deployed quickly, repeatedly, and flexibly. A 20-minute flight is sufficient to map a 10-ha area, and with advanced software, data processing can be completed in a few hours. Cameras designed for sUAS vegetation monitoring capture spectral information in visible and near infrared (NIR) bands, from which critical information on plant vigor and growth can be extracted. Users can easily customize the spatiotemporal resolution by adjusting the flight height and frequency. High-resolution digital surface models (DSMs) generated from sUAS data allow accounting for topo-edaphic conditions in forage production modeling and related analysis. A few pilot studies have explored the applications of sUAS technology in monitoring agricultural production for soybean [37], rice and wheat [38,39], barley [40,41], and mango [42]. However, there is currently limited research concerning the efficacy of sUASs for rangeland forage production quantification and modeling.
The primary aim of this paper was to develop LUE-based models to map forage production at very high spatial (30 cm) and temporal (daily) resolutions using multispectral data collected from sUAS, augmented with more frequent satellite data from PlanetScope satellites at 3-m resolution [43]. Specifically, our objectives included (1) imaging a Mediterranean-type semi-arid annual rangeland using sUAS and fusing the sUAS data with PlanetScope data to obtain daily sUAS resolution images; (2) investigating the connection between ground-based plant biophysical measurements and remote sensing vegetation indices, (3) building and evaluating forage production models to map daily rangeland production at 30-cm resolution, and (4) analyzing the predicted spatial and temporal patterns of estimated forage production to explore the relationship between forage production and its environmental and climatic drivers.

Study Site
Our study focused on a privately-owned, annual grassland (no trees or shrubs), located 56 km east of the coast in San Luis Obispo County, California (35.5N, 120.3W) (Figure 1), which has been grazed by beef cattle for more than 15 years. A 10-ha parcel was fenced in November 2016 to exclude grazing during the growing season. Following peak forage biomass growth, 50 cow-calf pairs grazed the fenced area for 30 days in the summer of 2017 to prepare the site for a repeated study during the 2018 growing season.
Soils were dominated by the loamy Balcom-Nacimiento complex [6]. The climate is Mediterranean with hot, dry summers and mild, moderately wet winters. Mean annual rainfall during 2001-2018 was 213 mm, mostly occurring from December through March. Annual rainfall at the site was 176, 139, 57, and 130 mm during the 2012-2016 drought years, respectively, and 287 and 123 mm during the study years of 2017 and 2018, respectively. Annual grasses and forbs at the site typically germinate in late fall with the onset of the rainy season (~November) and grow rapidly (March-April) to reach peak biomass

Field Measurements
In November 2016, 16 pairs of ECH2O 5TM sensors were installed to measure soil moisture and temperature across the study area ( Figure 1). Locations were selected based on a digital elevation model to capture the spatial variation of topography and soil conditions. At each location, duplicate sensors were deployed in the shallow rooting zone at a 7-cm depth (most roots occurred in the upper 15 cm soil layer). Data-loggers recorded sensor readings every 15 minutes. Three tipping bucket rain gauges recorded rainfall, and three time-lapse cameras followed forage germination and growth at the hilltop, south-facing slope, and north-facing slope positions ( Figure 1).
We measured forage production immediately after each of the 10 sUAS flight missions (except for 11/11/2016). Two replicate 30 cm × 30 cm quadrats at opposite cardinal angles were selected at a 1.5-m radius of the 16 soil sensors, resulting in two sample points per plot [44]. All aboveground plant biomass was harvested and dried at 60 • C for 48 hours before weighing to record dry biomass values. To further capture the spatial variation in topography and forage production, we classified the site into five categories based on clustering of topographical features including aspect, slope, flow accumulation, and elevation, using Iterative Self-Organizing Data Analysis Technique classification [45]. We selected an additional 20 plots for biomass clipping to ensure that each topographical cluster had at least eight plots. We collected 32~52 biomass clip-plot samples on each sUAS observation date and retained a total of 330 (200 and 130 in 2017 and 2018 growing seasons) biomass measurements after removing outliers (e.g., samples affected by strong rodent activity). For the cumulative NDVI & APAR and biomass analysis and LUE modeling, we further averaged the plots that had two samples near soil sensors, resulting 220 data points to work with.

sUAS Flights and Image Preprocessing
We integrated a MicaSense RedEdge camera with a 3DR Solo quadcopter for monthly aerial flights over the study area during the 2017 and 2018 growing seasons (November-April). MicaSense RedEdge is a multispectral camera with five spectral bands including blue, green, red, red-edge, and NIR. A sun irradiance sensor was included to measure band-specific incoming solar irradiance for radiometric correction. Flights were scheduled with the closest overpass dates of PlanetScope satellites when weather permitted it. A total of 10 missions were performed: six from 11 November 2016 to 30 April 2017, and four from 18 January 2018 to 14 April 2018.
We followed the same flight plan with a side-and front-overlap of 85% for all missions. Each mission was flown at 91-m (350-ft) above ground level (from the launching point) to acquire imagery at a 7.57 .9 cm resolution and cover the whole study area in less than 30 minutes due to the payload constraint. Flying speed was set at 7 m/s along a fixed direction parallel to the east/west site boundary ( Figure S1). Flights were conducted around solar noon to minimize the impacts of sun-angle variation throughout the season. Images of a calibrated white reflectance panel were recorded before and after each flight to calibrate raw pixel values to absolute reflectance values. To ensure the accuracy of the sUAS image geo-registration, we established eight permanent landmarks as ground control points (GCPs) evenly distributed at the corners (and center) of the study area. Coordinates of GCPs, forage clip-plot samples and soil sensors were recorded using a Trimble Geo 7x Series Handheld Data Collector.
Raw sUAS images were stitched and processed in Pix4dMapper Pro to generate orthomosaic maps of surface reflectance and DSM. The software converts the raw data in digital number (DN) to surface reflectance using images of the calibration panel, incoming sunlight irradiance, and camera parameters (e.g., ISO, aperture, shutter speed, and vignetting) recorded in the EXIF metadata [46]. DSMs were based on a dense point cloud generated from tie points that were automatically identified by the software. Mean average error (MAE) between GPS measurements and derived elevation was 3.5 cm across the study area.
A two-step geo-registration was performed on the time series of sUAS images to ensure spatial alignment of all aerial images. Images were first geo-registered in Pix4dMapper Pro using the eight GCPs to minimize the spatial mismatch when relating images to ground measurements. A mean Root Mean Square error (RMSE) of less than 3 cm was achieved. Reliability of this GCP-based geo-registration is highly dependent on the accuracy of the GPS. Since we did not upgrade our GPS to centimeter accuracy until the summer of 2017, we then used the Image to Image Registration function in ENVI 5.3 (using the 4/6/2017 image as the base image and second-order cubic) to align the images acquired in the 2017 growing season on different dates. In the Image to Images registration, we selected GCPs visually for each pair of images. The resulting RMSE was maintained at~3 cm, less than 0.3 pixels at the original 8-cm resolution for all pairs of images. The geo-registered 8-cm resolution surface reflectance maps were then aggregated to 30-cm resolution to match the plant biomass field measurement plots and used for further analysis, which resulted in a <10% geolocation displacement in the final 30-cm maps.
Standard reflectance calibration assumes similar incoming solar radiation for all pixels. However, the complex topography of the study site led to significant variation in illumination conditions among pixels at different topographic positions. We applied the C model [47] to correct pixel reflectance based on its illumination condition (IC) for each spectral band (supplemental online material, SOM). The C model has been applied for correcting Landsat satellite data [47][48][49][50], but its feasibility and accuracy in correcting sUAS has not been fully evaluated. We used our 30-cm resolution DSM derived from the sUAS data to perform the illumination correction. The illumination-corrected red and NIR reflectance values were then used to calculate sUAS NDVI.

High Resolution Satellite Imagery and Data Fusion
Relatively high temporal frequency of remote sensing imagery is needed to capture the rapid growth cycle of annual plants. We downloaded a total of 174 cloud-free PlanetScope (PS) orthorectified scenes (level 3A) from Planet at 3.125-m spatial resolution for the 2017 and 2018 growing seasons. With a constellation of around 120 CubeSats, PS aims to provide daily images of three visible bands and one NIR band for any place on Earth. We converted the DN to the Top of Atmosphere (TOA) reflectance, using the TOA reflectance coefficients embedded in the metadata XML for each individual scene. The NDVI was then calculated with red and NIR reflectance.
To compensate for the limited temporal acquisition of sUAS data, we combined PS satellite data with monthly sUAS data to obtain a complete time series of daily NDVI at the 30-cm sUAS resolution. Our goal was to predict a sUAS-resolution NDVI map U(x,y,t p ), for any particular date (t p ), using (1) two pairs of coincident sUAS and PS images during the nearest two sUAS flight dates on t b1 and t b2 and (2) a PS-based image on the prediction date t p . We followed the basic concept of the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [51] to take advantage of the complementary spatial and temporal information from PS and sUAS. We first took advantage of the PS's high temporal frequency and linearly interpolated temporally the available cloud-free PS NDVI images to obtain a continuous daily time series. To facilitate the pixel-based process, daily PS data were reprojected and bi-linearly resampled, denoted as P(x,y,t) to match the resolution, projection, and extent of the sUAS data. Since the study area has a relatively homogeneous landcover, we assumed that the temporal change from the daily PS NDVI images was similar at the 30-cm resolution: The sUAS-resolution NDVI for a given location (x,y) and day (t p ) was then estimated as a weighted sum of the two closest pairs of sUAS base NDVI images at 30-cm resolution ( Figure 2), each adjusted for temporal change derived from the PS NDVI time series: where t b1 and t b2 are the dates for the base image pair, and w 1 and w 2 are the temporal weights representing the contribution of each base image pair to the estimated NDVI, based on a linear function of the correlations between the base PS NDVI images on t b1 and t b2 and the PS NDVI image on t p (SOM).
where tb1 and tb2 are the dates for the base image pair, and w1 and w2 are the temporal weights

250
We tested the performance of our data fusion method using a leave-one-out method. Within We tested the performance of our data fusion method using a leave-one-out method. Within each iteration, we left out one sUAS NDVI image and predicted that image using the nearest neighboring image pairs and the PS NDVI image taken on the same day. The accuracy of the method was evaluated by comparing the predicted sUAS NDVI image to the observed sUAS NDVI image.

Empirical Statistical Analysis
We first examined the relationship between the measured biomass and two types of remote sensing metrics derived from the fused data, namely NDVI and APAR. NDVI was used as a proxy for net primary production and biomass [52][53][54][55]. Specifically, we used univariate linear regression (in R) to quantify the spatial variance in biomass across various plots explained by the corresponding coincident NDVI and the spatial and temporal variance explained when pooling data from all sites and dates together. When analyzing the coincident NDVI-biomass relationship, we did not average the paired samples taken close to the soil sensors, considering the limited number of samples on a single observation day. We investigated the relationship between biomass and cumulative NDVI integrated from the beginning of the growing season to the date of observation. Critical phenological days, such as the beginning (germination date) and end of the growing season were identified by fitting a logistic function (SOM) to the fused NDVI time series [56]. A similar statistical analysis was performed for cumulative APAR derived from the fused daily data and the California Irrigation Management Information System (CIMIS) spatial solar radiation.

Light Use Efficiency (LUE) Models
We estimated forage production based on LUE theory, as the cumulative product of APAR and LUE [22]: where, t p is the date of prediction, and t 0 is the corresponding germination date. APAR was estimated as the product of fPAR and PAR. We calculated fPAR as a function of NDVI [57]: where fPAR max = 0.95, fPAR min = 0.001, and NDVI max and NDVI min are the NDVI values corresponding to the upper and lower 2% of the NDVI histogram. NDVI min and NDVI max were set to 0.23 and 0.8, respectively, based on the NDVI distribution of sUAS data. PAR was calculated as 50% of the daily shortwave incoming solar radiation, available from the CIMIS Spatial dataset at a 2-km resolution [58], assuming a constant ratio of 0.5 for PAR over incoming solar radiation [59][60][61][62][63][64]. The CIMIS Spatial incoming solar radiation showed an R 2 of 0.99 when validated with CIMIS measured incoming solar radiation. The 2-km PAR was further sharpened to 30-cm resolution using the sUAS DSM to account for topography-induced illumination, and thus, PAR variations. We first generated daily solar radiation maps of the study area using the Area Solar Radiation function in the Arcpy package in Python, based on the sUAS-derived 30-cm DSM, latitude, and day of year [65]. Similar maps were derived for a simulated flat surface assuming DSM at the average elevation of the study area (491 m). For each 30-cm pixel and each day, the ratio of the daily 30-cm solar radiation from the actual DSM over the simulated flat DSM was used to multiply the corresponding CIMIS spatial solar radiation to derive solar radiation and PAR at 30-cm resolution.
We developed two semi-empirical statistical models to estimate LUE, depending on availability of the input data. All parameters in LUE models were optimized with field measurements. We first parameterized LUE as a function of topographic variables (Equation (4), Model I), which were derived from the DSM generated from sUAS data: LUE 0 is the LUE at the optimal topographic position. LUE 0 is estimated by scaling f(Topo,x,y) from 0 to 1. We used topographic variables to drive an empirical scalar because we did not have the spatial soil moisture and air temperature data that are traditionally used for LUE down-regulators. We explored a suite of topographic features including cosine aspect, slope, flow accumulation, IC, and elevation, and enumerated the possible combinations of these variables in the equation.
We further added time-varying controls of soil moisture (W) and air temperature (T) on LUE, in addition to significant topographic features (Equation (5), Model II), as shown below: where g W, x, y, t p , p T, x, y, t p , and h(Topo, x, y) range from 0 to 1. Forms of the W and T scalar functions were set to the same as those in Reference [22] but coefficients were optimized using field-measured biomass. The temperature (T) stress scalar was based on deviation from optimal air temperature (T opt ): where a is the coefficient to be calibrated, T a is the air temperature ( • C), and (x, y, t p ) denotes the location and time. T opt is here defined as the average temperature during the month with the highest NDVI. We used daily mean air temperature from the 2-km CIMIS spatial dataset. The T scalar decreases from 1 to 0 as temperature deviates from T opt , with the rate of decrease represented by a.
The water stress (W) scalar represents the down regulation of light use efficiency under drought conditions, and was parameterized as a function of soil moisture (M) (Equation (7)): where b is the coefficient to be calibrated; M is the derived soil moisture, represented by volumetric water content (VWC); and WP is the wilting point. W ranges from 0.5 to 1 as soil moisture varies from the wilting point to field capacity. The wilting point was set to a water content of 0.14 cm 3 cm −3 determined by averaging soil moisture content in June during the dry season when grass and forbs had desiccated and soil moisture was at its water-year minimum [66]. For operational mapping purposes, we required soil moisture for every pixel. We estimated the daily soil moisture using a simple bucket model (SOM) following Reference [22], based primarily on precipitation data from PRISM [67] and potential evapotranspiration from CIMIS Spatial [58]. The topographic control term is a linear function of elevation and IC for Model II. Elevation and IC were selected because they had the highest statistical significance for minimizing differences between measured and predicted biomass among topographic variables. As a comparison to Model II, we also built a Model III by replacing the air temperature and derived soil moisture with measured soil temperature and soil moisture in the LUE parameterization. As for the topographic variable, we selected elevation because it had the highest statistical significance on minimizing model error.
We optimized the model coefficients using the Stochastic Gradient Descent (SGD) method [68,69]. We used 70% of the data for training and the remaining 30% for validation of the SGD LUE optimization. Model performance was evaluated by comparing model estimates with measured biomass values in the validation dataset. We used R-squared (R 2 ) and RMSE to quantify the uncertainty of models for estimating forage production. The LUE scalars and their corresponding coefficients may be site-specific as they were specifically optimized for our study area.

Forage Production Mapping and Patterns
We implemented Model II to generate daily forage production maps for the study area. The maps were compared against the concurrent RGB images from sUAS to determine if prediction captured a similar spatial pattern. To further study the interaction between topography (aspect and slope) and forage production, we calculated the predicted biomass for topographic zones divided by aspect (north and south) and slope (flat, moderate, and steep) and performed zonal statistical analysis.

Terrain Correction
The wavelength-specific C model significantly removed the topographic varying illumination effect present in the original sUAS images ( Figure 3). For example, pixels on the north-facing slopes, when displayed in true color composite, appeared very dark due to terrain shadows, especially for early season images when the sun angle was relatively low, e.g., 16 January 2017. This shadowing effect was minimized after applying the correction. Raw reflectance was significantly correlated with IC (p < 0.01 in most cases), but the correlation was minimized and not significant after correction ( Figure 4). For example, R 2 was reduced from 0.33 (p < 0.01) to 0.12 (p = 0.02) in the red and from 0.80 (p < 0.01) to 0.01 (p = 0.5) in the NIR for January 2017 imagery ( Figure 4). The topographic illumination effects were found to be more pronounced in the NIR band and early in the growing season.

sUAS and PlanetScope Data Fusion
We combined all available 30-cm sUAS NDVI images and satellite NDVI images from the two growing seasons to generate daily 30-cm NDVI time series from 11 November 2016 to 30 April 2017 and from 1 January to 14 April 2018. The fused daily NDVI curve had a similar magnitude to those from sUAS imagery, and a quantitative comparison showed that the predicted 30-cm NDVI agreed well with the original sUAS NDVI images (R 2 = 0.87 and RMSE = 0.06) (SOM, Section 2). The fused time series followed similar temporal patterns with those from PS data. For a selected plant pixel within a relatively homogeneous patch, Figure 5 shows three sets of NDVI values from PS and sUAS sensors at 30-cm resolution and data fusion. The sUAS NDVI values were typically higher than those from the concurrent PS imagery, e.g., 0.65 vs. 0.46 on 6 April 2017. This difference was most likely due to lower atmospheric effects on surface reflectance from the low-altitude sUAS flights. The fused time series was able to conserve the sUAS NDVI values, e.g., 0.65 on the same day. Due to its higher temporal frequency, the PS NDVI time series captured more detailed temporal dynamics of plant growth than the less frequent sUAS snapshots. Both fused and original PS NDVI time series in 2018 showed a unique bimodal shape, with peaks on January 25 th and March 19 th . This variation was consistent with a gap in precipitation that greatly reduced plant growth in February, but growth resumed following late March rainfall. The smoothed and fused NDVI time series successfully captured the rapid changes in phenology, while preserving the magnitude of NDVI values with minimal atmospheric contamination ( Figure 5).

Relationships Between Remote Sensing Metrics and Biomass Measurements
Daily fused NDVI was positively related (with varying significance) to biomass measurements (mean = 1122 kg/ha; range = 17 -4483 kg/ha) taken on the same day over all the clipping sites, with R 2 ranging from 0.01 to 0.33 on various sUAS image acquisition dates ( Figure S4). The highest R 2 , 0.33 (p < 0.05, N = 21), was found for February 2018. The varying R 2 among different dates suggests limited performance for using instantaneous NDVI to predict biomass over the growing season. NDVI values were typically lower in 2018 than those in 2017 for the same month, consistent with differences in corresponding biomass measurements.
When integrated from the germination date to the date of observation (using the fused time series), the cumulative NDVI explained 68% of the variance in measured biomass among all field sampling plots during the two growing seasons (Figure 7a), and the cumulative APAR explained 69% of spatial and temporal variance (Figure 7b) (N = 220). The linear regression model, based on 70% of the randomly selected data points, showed that cumulative APAR predicted biomass with an R2 of 0.67 ± 0.06 and a RMSE of 631 ± 82 kg/ha, when compared with the remaining 30% of the data (N = 66). Similarly, predictions from cumulative NDVI showed an R2 of 0.67 ± 0.05 and an RMSE of 627 ± 73 kg/ha. However, the empirical NDVI-only and APAR-only methods may have large uncertainties when extrapolated to other areas or time periods.

Forage Production Models
Optimization of the topography-based LUE model (Model I) with 70% of randomly selected data resulted in the following forage production model: The soil moisture scalar showed that decreasing moisture down-regulated LUE. The W scalar increases logistically from 0.5 to 1 as soil moisture changes from wilting point (θ v = 0.14 cm 3 cm −3 ) to field capacity (θ v = 0.29 cm 3 cm −3 ). The T scalar dropped to 0.5 when air temperature was 1 • C below or above the optimal temperature of 10.33 • C. This full model resulted in improved accuracy for predicting biomass with greater explanation of variance (81%) and lower RMSE (542 kg/ha) compared with the topography-based LUE model (Figure 8b). This suggests that time-varying controllers of soil moisture and air temperature captured environmental stress impacts on forage production.
When using field-measured soil moisture and soil temperature, Model III: predicted biomass with an RMSE of 472 kg/ha and R 2 of 0.77 ( Figure S5).

Forage Production Mapping and Patterns
We summarized forage production across the study area based on daily biomass maps generated from Model II. Forage production increased gradually during the growing season, but followed an unstable trajectory with bursts and plateaus (Figure 9). Most rapid growth was from 20 March to 2 April 2017 (averaging 107 kg/ha/d) and from 10 March to 19 March 2018 (averaging 15 kg/ha/d). The 2017 growing season had higher overall growth rates than 2018, which resulted in a mean peak standing biomass of 3216 (±678SD) kg/ha versus 1054 (±374SD) kg/ha, respectively. This large difference in peak standing biomass was linked to higher precipitation in 2017 (287 mm) compared to 2018 (123 mm). The standard deviation of biomass for the sUAS flight dates (Figure 9) showed increasing spatial variability in biomass as the growing season progressed.

451
We summarized forage production across the study area based on daily biomass maps 452 generated from Model II. Forage production increased gradually during the growing season, but 453 followed an unstable trajectory with bursts and plateaus ( Figure 9). Most rapid growth was from mean peak standing biomass of 3216 (±678SD) kg/ha versus 1054 (±374SD) kg/ha, respectively. This  Predicted biomass maps at 30-cm resolution provided a visual representation of spatial variation across the landscape (Figure 10). At peak biomass in 2017 (Figure 10d), high biomass production was predicted at low elevation ( Figure 10i) and high IC (Figure 10j) regions, whereas in 2018, high biomass production was predicted in low elevation and low IC regions. Lower elevation regions often have higher soil moisture because lower slope positions receive runoff and throughflow from upper hillslope positions. IC values are positively related to PAR with higher PAR leading to higher biomass when plants are not water stressed. However, when water is limited, high PAR could contribute to lower biomass production (i.e., greater plant water stress). The maps also captured human disturbance on biomass production, such as the road having very low biomass compared to the surrounding area.

487
In this pioneering study, we developed two remote sensing data-driven LUE models that 488 accurately predicted daily forage production at centimeter-scale resolution. We also streamlined 489 preprocessing and fusing approaches for raw sUAS data with high-resolution satellite remote sensing 490 data to obtain high-spatiotemporal data. We successfully implemented preprocessing and fusing forage production maps, we observed unique patterns in productivity related to different

498
This study demonstrated a successful application of the C model for correcting sUAS data for 499 illumination effects. Reflectance in visible bands was generally less affected by terrain than NIR and 500 red-edge bands. The same phenomenon was found for topographic illumination effects on Landsat 501 images [50]. The authors concluded that the lower correlation for visible bands was caused by their 502 lower reflectance rate and a more significant atmospheric scattering effect on visible bands than 503 longer-wavelength bands. From the temporal analysis, the topographic illumination effect intensity 504 was not consistent over the entire growing season for a given band, but rather diminished as the 505 growing season progressed. In addition to changes in vegetation density, the diminishing trend is 506 probably related to changes in sun angle. Since the study site is at 35.5N, the solar-noon sun elevation 507 angle increases from 31˚ in December to >80º in June. Therefore, the theoretical intensity of 508 illumination effects should decrease as the growing season progresses.

509
Notably, the R 2 between red band surface reflectance and IC was very low (0.02) for the March 510 flight, and the corresponding corrected R 2 (0.06) was higher than the original R 2 . Although both R 2 511 values were small, this raises caution that applying the C model may lead to overcorrection when the 512 illumination effect is weak. Therefore, a test of the illumination effect intensity, for example, a 513 correlation analysis on the surface reflectance and IC, is recommended before applying illumination 514 correction to models. In this study, since the corrected red band surface reflectance-IC R 2 value was 515 very low, we used the corrected red band reflectance to maintain temporal consistency among data.

Discussion
In this pioneering study, we developed two remote sensing data-driven LUE models that accurately predicted daily forage production at centimeter-scale resolution. We also streamlined preprocessing and fusing approaches for raw sUAS data with high-resolution satellite remote sensing data to obtain high-spatiotemporal data. We successfully implemented preprocessing and fusing methods and frameworks (e.g., the C model for illumination correction and the simplified STARFM) that were developed for satellite remote sensing data on the sUAS data. From the model predicted forage production maps, we observed unique patterns in productivity related to different environmental and climatic drivers. The following discussion section expands on additional findings in sUAS data preprocessing (Section 4.1) and fusing with PlanetScope data (Section 4.2), and the response of forage plants to precipitation (Section 4.3) and soil moisture and temperature (Section 4.4).

Variations in the Illumination Effect
This study demonstrated a successful application of the C model for correcting sUAS data for illumination effects. Reflectance in visible bands was generally less affected by terrain than NIR and red-edge bands. The same phenomenon was found for topographic illumination effects on Landsat images [50]. The authors concluded that the lower correlation for visible bands was caused by their lower reflectance rate and a more significant atmospheric scattering effect on visible bands than longer-wavelength bands. From the temporal analysis, the topographic illumination effect intensity was not consistent over the entire growing season for a given band, but rather diminished as the growing season progressed. In addition to changes in vegetation density, the diminishing trend is probably related to changes in sun angle. Since the study site is at 35.5N, the solar-noon sun elevation angle increases from 31 • in December to >80 • in June. Therefore, the theoretical intensity of illumination effects should decrease as the growing season progresses.
Notably, the R 2 between red band surface reflectance and IC was very low (0.02) for the March flight, and the corresponding corrected R 2 (0.06) was higher than the original R 2 . Although both R 2 values were small, this raises caution that applying the C model may lead to overcorrection when the illumination effect is weak. Therefore, a test of the illumination effect intensity, for example, a correlation analysis on the surface reflectance and IC, is recommended before applying illumination correction to models. In this study, since the corrected red band surface reflectance-IC R 2 value was very low, we used the corrected red band reflectance to maintain temporal consistency among data.

Fusing Satellite and sUAS Data
We found that the PS NDVI values were generally lower than those from sUAS, especially when NDVI was higher than 0.5. This results from PS NDVI being calculated from the Top of Atmosphere (TOA) reflectance while the sUAS NDVI is determined from the surface reflectance. Several studies document a lower TOA NDVI compared to NDVI at the surface [70,71] due to the atmospheric scattering effect. In the simplified STARFM, we used the temporal change of PS TOA NDVI, which reduced the noise from TOA reflectance, while the absolute value relies heavily on the sUAS. Smaller differences between December and February sUAS instantaneous NDVI and PS NDVI datasets may result from the use of a sun light irradiance (SLR) sensor on sUAS flights after March 2017. We initially relied on a white calibration panel for radiometric and atmospheric correction when preprocessing data from the first two flights before installation of the SLR sensor.
Data fusion of sUAS and high spatial-resolution images has great potential to lower the cost of operating frequent sUAS flights, especially given increasing availability of high-spatial resolution data (e.g., RapidEye data and PS data by Planet and WorldView data by DigitalGlobe). However, much of the high-spatial resolution satellite (including PS) data is only available in DN or TOA due to the lack of shortwave NIR bands for atmospheric correction. Therefore, it is important to develop robust data fusion methods that work on multi-source data with systematic differences. Future release of atmospherically-corrected PS surface reflectance products, i.e., using the concurrent MODIS data, will further reduce uncertainties in sUAS-satellite data fusion methods. Moreover, a more sophisticated fusion method is needed to take into account the spatial details embedded in sUAS data for the temporal interpolation of PS data.

LUE Parameterization
We observed improvements in model accuracy when including additional model terms. RMSEs of forage production estimates from APAR-only, Model I, Model II, and Model III were 624 kg/ha, 567 kg/ha, 542 kg/ha, and 472 kg/ha, respectively. We selected Model II to map forage production because it had the highest accuracy and does not require ground measured soil moisture and temperature. Model III, which uses locally measured soil moisture and soil temperature, achieved a higher accuracy than Model II, which uses estimated soil moisture and coarse-resolution air temperature. We believe that the main source of error in Model II is from the use of a simple bucket model for estimating soil moisture and/or the coarse resolution of air temperature.
Comparing the stressing scalars in Model II and III, the W scalar functions have very similar controls on LUE; however, the T scalar behaves quite differently. In Model II, the T scalar is driven by air temperature and is more sensitive to deviations from optimal air temperature. The estimated LUE 0 across the three models showed large discrepancies, ranging from 0.32 to 1.8 g/MJ APAR with Model I having the highest LUE 0 and Model II the lowest. Over the years, maximal LUE has been estimated using different methods for different biomes [72]. However, maximal LUE used in existing LUE models varies a lot [31] because researchers used different methods to estimate this metric [72]. In addition, maximal LUE also varies across scales. Currently, the number of LUE studies at the watershed scale is very limited. We expect to see increasingly more relevant studies on LUE modeling with advances in sUAS technology.
Using measured soil temperature and moisture in Model III achieved the highest accuracy among the three models, but this model cannot be applied to directly generate biomass maps because it requires point measurement data to parameterize input variables. As of now, daily soil moisture maps at high spatial resolution are not available, but with the continuous expanding constellations of satellites, new remote sensing products are being produced that could provide daily spatial patterns for soil moisture. For example, a past study demonstrated the potential for retrieving soil moisture from C-band synthetic aperture radar (SAR) data [73]. Additional studies have developed algorithms for mapping high-resolution (<1 km) soil moisture locally using the sentinel-1 C-band SAR data at a 10-m resolution [74,75]. Therefore, inclusion of near real-time soil moisture estimates may be forthcoming as new remote sensing technologies emerge.

Response of Forage Production and Plant Phenology to Moisture
Moisture is the primary controlling factor of forage production in Mediterranean annual range systems [6,31,76]. The different precipitation regimes in the 2017 and 2018 growing seasons triggered very different biomass-aspect relationships. South-facing slopes were expected to have a lower peak biomass than north-facing slopes because of increased temperatures leading to higher ET and lower soil moisture on south-facing slopes. However, we observed a higher peak biomass on south-facing slopes in the 2017 growing season (Figure 10b). This resulted from the high precipitation (287 mm) that was well distributed throughout the 2017 growing season. With sufficient precipitation to maintain soil moisture content ( Figure S6), plants experienced little water deficit. With sufficient water supply, PAR became the primary limiting factor. Therefore, south-facing slopes, where more solar radiation is received, were able to produce higher peak biomass than north-facing slopes. In contrast, the 2018 growing season received only 108 mm of precipitation with large gaps between rainfall events. The majority of the precipitation was received in January, March and early April resulting in an extremely dry February (Figures S5 and S6). This precipitation distribution hindered plant growth on the radiation-rich south-facing slopes by intensifying the soil moisture deficit. Therefore, a lower peak biomass was observed on south-facing slopes in the 2018 growing season.
Time series biomass maps provided insights on plant phenology changes in response to different precipitation regimes during the 2017 and 2018 growing seasons. Germination and growth can be very different interannually in annual range systems. Peak growth is highly dependent on the amount and timing of precipitation. The November 2017 rainfall resulted in germination (11/11/2016) of the plants about 2 months earlier than the 2018 growing season (1/1/2018) (Figures 5 and 12). The predicted germination dates are close to actual gemination dates (11/6/2016 and 1/13/2018) captured by the time lapse camera. Although the two growing seasons had a 1.5-month difference in their germination date, they both had peak growth in April. We predicted 4/6/2017 and 4/14/2018 as the peak growth dates for the 2017 and 2018 growing season, which were close to the camera observed peak growth dates (4/12/2017 and 4/25/2018). In addition to having a longer growing season, plants also accumulated biomass faster in 2017 than in 2018 ( Figure 9). This higher biomass accumulation rate is largely due to the higher and more evenly distributed precipitation pattern during 2017. Plant growth on south-facing slopes was very sensitive to soil moisture availability. During the 2017 growing season and in January 2018, when moisture was sufficient, plants on south-facing slopes accumulated biomass faster than those on north-facing slopes. However, when soil moisture became limiting, for example, starting in February 2018, biomass accumulation on south-facing slopes was distinctly lower than north-facing slopes. water supply, PAR became the primary limiting factor. Therefore, south-facing slopes, where more 576 solar radiation is received, were able to produce higher peak biomass than north-facing slopes. In

Conclusions
This pilot study demonstrates the synergistic use of complementary sUAS and satellite remote sensing technology to map daily forage production at 30-cm resolution on a delayed grazing rangeland. Remote sensing-based APAR estimates were the primary driver for both LUE-based models developed here. LUE was optimized as a function of elevation and illumination conditions in Model I, achieving an R 2 of 0.70 and an RMSE of 567 kg/ha when comparing predicted forage production with measured forage biomass. By further incorporating moisture and temperature stress terms to the LUE parameterization (Model II), Model II predicted forage production with higher accuracy (R 2 = 0.81) and lower RMSE (542 kg/ha). Both LUE approaches showed improved prediction accuracy compared to the univariate linear regression with cumulative APAR. Finally, the inclusion of measured soil moisture and soil temperature into the LUE model (Model III) achieved the lowest RMSE (472 kg/ha) with a similar R 2 of 0.77.
The fusion of high spatial resolution sUAS data and high temporal PlanetScope satellite imagery enabled us to generate daily forage production maps at a 30-cm resolution. The maps allowed the analysis of landscape-derived zonal statistics to assess (a) interactions between topography and biomass and (b) the response of forage production and plant phenology to changing precipitation regimes. Our analysis showed higher forage production and biomass accumulation rates in landscape positions experiencing less environmental stress (e.g., soil water deficit), with several differences occurring between wet and dry years. The forage production maps can be implemented in decision support tools to help ranchers better anticipate weather-driven changes in forage production and optimize their decisions on proactive practices, such as stocking conservatively, resting pasture, and incorporating yearling cattle. This study provides basic tools that can be further developed and scaled to statewide regions to provide near real time forage availability delivered to users via internet apps.