Comparing Two Methods of Leaf Area Index Estimation for Rice ( Oryza sativa L.) Using In-Field Spectroradiometric Measurements and Multispectral Satellite Images

: This work presents a remote sensing application to estimate the leaf area index (LAI) in two rice ( Oryza sativa L.) varieties (IDIAP 52-05 and IDIAP FL 137-11), as a proxy for crop performance. In-ﬁeld, homogeneous spectroradiometric measurements (350–1050 nm) were carried in two campaigns (June–November 2017 and January–March 2018), on a private farm, TESKO, located in Juan Hombrón, Coclé Province, Panama. The spectral ﬁngerprint of IDIAP 52-05 plants was collected in four dates (47, 67, 82 and 116 days after sowing), according to known phenological stages of rice plant growth. Moreover, true LAI or green leaf area was measured from representative plants and compared to LAI calculated from normalized PlanetScope multi-spectral satellite images (selected according to dates close to the in-ﬁeld collection). Two distinct estimation models were used to establish the relationships of measured LAI and two vegetational spectral indices (NDVI and MTVI2). The results show that the MTVI2 based model has a slightly higher predictive ability of true LAI ( R 2 = 0.92, RMSE = 2.20), than the NDVI model. Furthermore, the satellite images collected were corrected and satellite LAI was contrasted with true LAI, achieving in average 18% for Model 2 for MTVI2, with the NDVI (Model 1) corrected model having a smaller error around 13%. This work provides an important advance in precision agriculture, speciﬁcally in the monitoring of total crop growth via LAI for rice crops in the Republic of Panama.


Introduction
Rice (Oryza sativa L.) is an essential grain in the primary diet of many people worldwide, especially for inhabitants of developing countries [1]. About 23% of the calories consumed by the world's population come from rice [2]. Rice is mostly produced and consumed within the Asian continent. This continent is also responsible for more than 90% of its world production [2]. In Panama, rice was declared an essential crop for food security, through Law 17 of 22 February 2018 , as this is the main product of the basic food basket in the country (Official Gazette No. 28471-B). studied the effect of panicle presence on LAI estimation for monitoring rice cultivation at the harvest stage.
This article analyzes rice cultivation's spectral signature to obtain the bands that best correlate with the LAI and uses two distinct models for its prediction. The purpose of this work is to contribute to knowledge regarding the use of spectral data from satellite images for the estimation of biophysical variables through remote sensing; specifically, it seeks to establish the basis for the development of models that incorporate vegetation indices for the estimation of the leaf area index, which is applied in the description of the state and productivity of rice cultivation, using PlanetScope images.

Characterization of the Study Area
The study site is located on the TESKO S.A. farm, in Juan Hombrónn, Coclé Province, at coordinates 8 • 19 N, 80 • 13 W, approximately 150 km west of Panama City (see Figure 1 for detailed location). The farm is located in Juan Hombrón, district of Antón, province of Coclé. This location is part of the Dry Arc region (the Central Pacific Hydrologic Region of Panama). This region has an average precipitation of 1400 mm per year, in which dry months precipitation could be as low 10 mm/month and in wetmonths up to 500 mm/month [31]. According to Köppen, this region is considered to have a dry-humid tropical climate or tropical Savannah climate. Juan Hombrón's has an economy based on agricultural activities, with agri-businesses focusing on crops such as: Sugar cane, coffee, corn, rice, beans. In addition, fishing, shrimp farming and livestock are grown in this area.

Cultivar Selection
Originally, IDIAP 52-05 variety was chosen for all experiments in this study due to having a good milling performance and a superior industrial quality. IDIAP 52-05 fulfills its life cycle between 116-124 days after sowing. Plants reach maximum heights of between 85-110 cm and are adapted to both irrigation and rain-fed conditions. More importantly, it is known to be resistant to Pyricularia sp. among other pathogens [32].
Therefore, the base unit of analysis for this study were crops of the rice variety IDIAP 52-05 grown on the TESKO S.A. farm. However, it should be stated that management decisions made by the owners of the farm and therefore out out of our hands, impacted the course of this research project and forced the group to use IDIAP 137-11, a similar rice variety to IDIAP 52-05, in order to validate the results. The IDIAP FL 137-11 variety has an early-intermediate vegetative cycle, with 107 to 129 days from sowing under rainfed conditions and 108 to 114 days under irrigated conditions. It has a height that oscillates between 83 and 114 cm [33]. This is of important to establish, since IDIAP FL 137-11 crop with 77 days after sowing, was used for validation of the models.

Sampling Dates
Rice crops are often divided by the stage in their growth phase. In literature, rice phenology is divided into three main phenological phases (as described in [33,34]), as follows: • Phase #1: vegetative, characterized by starting with germination and concluding with panicle initiation (days 1 to 55 after sowing). • Phase #2: reproductive, characterized by panicle initiation and concluding in flowering (days 55 to 105 after sowing). • Phase #3: maturity or ripening, which is characterized by grain filling and extending until maturity (days 105 to 120 after sowing).
As the rice variety under study develops for about 4 months. The sampling dates, were observed to line up with phenological stages, as follows: (1) vegetative, with a sampling date 47 days after sowing, (2) reproductive with with two sampling dates 67 and 82 days after sowing and (3) maturity with a sampling date 116 days after sowing. This linking is necessary in order to explain their evolution according to management operations for rice crops.
It should be noted that spectral signature at the beginning of the crop on the early vegetative stage (13 days after sowing) were collected. However, as this spectral signatures suffer from suffering from having a mayor reflectance component coming from the soil (due to small leaf area). Measurements of this type are described in length in Sanchez-Galan et al. [35,36]).

In-Field Spectral Signature Collection
For the capture of spectral signature data, preferential sampling was chosen. That is, taking several replicas of the spectral signature of the coverage class previously identified at homogeneous points, which characterize and represent that class. A spectroradiometer (GER 1500 from Spectra Vista Corporation) covering a range of 350-1050 nm was used, with a field of view of 8 • . Before each signature shot, a reference measurement is captured on a Spectralon (a calibrated white semi-Lambertian surface). Figure 2 visually represents the team acquiring a in-field with portable spectroradiometer. Spectroradiometric measurements were made between June and September 2017 and between January and April 2018. All measurements were acquired in an specific period of the day, either between 10:00 a.m. and 11:30 a.m. or between 1:00 p.m. and 2:30 p.m. This ensured controlled reflactance variability and signal quality. In all cases, a blank (reference measurement), were taken between samples rapidly as the sunlight and clarity conditions changed appreciably. Each signature measurement was georeferenced using a Garmin eTrex GPS. Finally, the data were downloaded as text files with ".asc" extension and saved in specific folders. File names were saved with a code indicating the rice variety and the measurement date.

Transformation of the On-Site Spectral Signature to a Satellite Spectral Signature
Since the equipment manufacturer adhered to the ISO 80000 standard https://www. iso.org/obp/ui/#iso:std:iso:80000:-7:ed-2:v1:en:en%00 (Last accessed on: 2 May 2023) and provided the measurement in units of W × cm −2 × nm −1 × sr −1 × 10 −10 , it was necessary to make a transformation of spectral signature given by the spectroradiometer to the ones given by the remote sensor from PlanetScope, operated by Planet Labs [37]. The solar irradiance was analyzed after a transformation to W × m −2 for each band used, using the Equations (1) and (2): (1) where SI is solar irradiance, BR is blank radiance TI is target irradiance, and TR is target radiance. λ n corresponds to the wavelengths in nm and m is the number of wavelengths that the spectroradiometer has within each of the PlanetScope satellite bands (ρ = {blue, green, red, and infrared}). Ω corresponds to the solid angle in steradians (sr), which is related to the field of view by the Equation (3): For the particular case the in-field equipment had an α = 8 • , corresponding to a solid angle of Ω = 0.0153 steradians. Finally, the 10 −6 factor results from the conversion cm 2 to m 2 . Thus, the reflectance of a band ρ is obtained by the relationship shown in Equation (4): Moreover, on-site spectral signals were re-calculated to take into account the Relative Spectral Response (RSR) filters, that is, the specific Planetscope satellite relative sensitivity to radiance at various parts of the electromagnetic spectrum. Filters were obtained directly from Planetscope website for the Dove Classic Satellite.
As the satellite images used in this study were acquired within the 2017-2018 period. Our images belong to the "Dove classic" instrument. Therefore, the sensitivity files were be obtained from a PlanetScope forum https://support.planet.com/hc/en-us/articles/360014 290293-Do-you-provide-Relative-Spectral-Response-Curves-RSRs-for-your-satellites-(Last accessed on: 2 May 2023), with the code 0fxx/Dove-Classic_10xx. The RSRs have a value between 0 to 1 and are unit-less, since they are described relative to the peak response; for out of this case, they are shown in Figure 3.
Equation (5), describes the relationship between the reflectance from the spectroradiometer and the reflectance for each the wavelength of each sensor, including the RSR filter.
where R band , is the reflectance at a range of wavelength (e.g., red band). R i , is the reflectance at a specific wavelength, and RSR i , it the relative response value at the specific wavelength. This equation have been successfully used in Melillo et al. [38], Trischenko et al. [39,40] and more importantly in Agapiou et al. for the same GER 1500 Spectra Vista Corp. instrument as used for our measurements [41].  Figure 4 shows the process of collection of spectral signatures. In addition, the sampling area selected for True LAI determination. Representative points of the spectral signature of the crop were selected within the plot. Five random samples were cut in triplicate at days 47, 67, 82, 97 and 116; and collected in an area of 25 × 25 cm 2 for the complete cultivation cycle (about 4 months). As demonstrated in Figure 4A, the sampling process consisted of cutting the plant's entire aerial sample with scissors and placing it in a closed bag. As shown in Figure 4B the sampling area was separated with a wooden stick frames and threads, for consistency and repeatability. Samples collected over the experimental area were taken to the laboratory to estimate the green leaf area (dried leaves were not included for the calculation of the area). For this study, this measurement is considered as true LAI. Then, each rice leaf was mounted on a flat white cardboard surface, taking into consideration a proportional separation between samples. A calibrated ruler was placed on the cardboard and used as a standard to have the profile object of known dimensions used later for comparison with the leaves. Finally, the mounted leaves were photographed from the top, and the image was saved accordingly for later estimation.

Measurements of Green Leaf Area
The photographs were analyzed with the ImageJ program version 1.51t by selecting a mask to obtain the borders of the leaf finally the shaded area corresponding to the leaf area was calculated [42]. This procedure was repeated twice for selected crop dates. Figure 5, provides a visual explanation of the measurement procedure. Figure 5. True LAI Determination via the ImageJ software. Leafs were arranged in the controlled background surface (A), then the reference and the green leafs were selected and highlighted area confirmed in the software interface (B). Finally, the green area was calculated (C).

Modeling the Relationship between In-Situ Reflectance and LAI
Regarding multispectral satellite imagery, PlanetScope images were selected, as previously agreed with the funding agency. Selected images had a 3 × 3 m resolution. The satellite bands used in this study were: ρ blue = 455-515 nm, ρ green = 500-590 nm, ρ red = 590-670 nm and ρ N IR = 780-860 nm., respectively. Using the image search engine of Planet Labs, a private company that processes and distributes Satellite Images PlanetScope, images were identified, selected, and downloaded according to the georeferenced points. Each image was actually a folder with three images. Images that in their name include the suffix "_SR" have radiometric and atmospheric corrections.
These georeferenced satellite images needed to match two criteria: (1) to be in dates as close as possible so that they corresponded to the age of rice cultivation (in days) and (2) those with the lowest percentage of cloud cover. If the captured image did not coincide with the date of measurement in the field or had a significant cloud percentage (as is typical for Panama), images of the nearest day in the same week for the age of the crop were downloaded. These caused a partial time discrepancy between field measurements and satellite acquisition date.
Regarding the images' quality, they all have atmospheric corrections to obtain the groundlevel reflectance (surface reflectance). All images used had two predetermined corrections made by the supplier: atmospheric correction algorithm: "6Sv2.1" and aerosol model: "continental".

Spectral Estimation of LAI
Two models were used to understand the relationship between true LAI and in-situ reflectance during the vegetative and reproductive phases. Model #1 is directly related to NDVI and Model #2 with the MTVI2. Both, Models #1 and Model #2, were first estimated using Equations (6)- (8), found in literature. Model 1 was first suggested by Nguyen et al. [43] and Model 2 by Haboundane et al. [13], but later estimated to fit the acquired data. where Model validation was made by adding, for all points considered, the percentual errors, by using Equation (9): where the true value is represented by the True LAI and the calculated value is the one estimated by each model.

Calibration Adjustment to PlanetScope Satellite Imagery
Additional radiometric calibration settings were employed relative to initial corrections and referenced to the irradiance and reflectance measured in the field (as a function of the day and time). In other words, the measured in-field spectra were the best approximation of the spectral signature with the slightest error due to the intervention of cloud reflection or infrared re-absorption. The correction factor is the ratio of the field reflectance to the satellite reflectance using the band ranges delimited by the PlanetScope sensors.
The correction factor was obtained with the data of four different days from Juan Hombrón plot: 1.
25 January 2018 (field and satellite measurement). For all four days the correction factor was estimated as the ratio given in Equation (10): Correction factor = field reflectance satellite reflectance (10) A correction factor regression with solar irradiance W m 2 per band was performed to determine the change in correction factor due to changes in solar irradiance. Finally, reflectance is calculated using Equation (11): Corrected reflectance = satellite reflectance × correction factor (11)

Software
All satellite imagery was processed with the QGIS 3.10 software Package [44]. The R program was used for the statistical analysis of data [45], including spectral band integration and figure creation. Statistical significance was established with a cut-off of 0.05 (p < 0.05). Second order non-linear polynomial regression needed for calibration betweeen in-field and satellite images were calculated using the nls https://cran.r-project. org/web/packages/nls2/index.html (Last accessed on: 2 May 2023) package in R. The Pseudo R 2 estimates for the non-linear regression [46] were calculated using the aomisc (https://rdrr.io/github/OnofriAndreaPG/aomisc/ (Last accessed on: 2 May 2023) package from R stastical software. Figure 7 shows the spectral signatures of the IDIAP 52-05 variety for four sampling dates. Figure 7A provides a representation of vegetative phase on day 47 after sowing. Figure 7B,C represents the reproductive phase on days 67 and 82 after sowing, respectively. Figure 7D, shows the maturity phase, specifically on day 116 after sowing. This figure shows the evolution of the spectral signatures over time and, more importantly, its spectral separability. One can notice that the first and last phases have a lower overall reflectance, contrasting with the reproductive having higher reflectance, with notable reflectance attributed to the plant's flowering.

True LAI Estimation
As it was previously described in Section 2.5, true LAI was calculated via image analysis of rice plant leaves. The leaves were cut out, and then the ImageJ program was used to count the green area of the leaves in the photograph scaled with an object of known dimension. Table 1, shows the estimated average true LAI for the IDIAP 52-05 variety at different days after sowing. The Mean Absolute Deviation (MAD) was calculated to assess the dispersion values in LAI measurements for each sowing date. For the calculation of the MAD Equation (12), was used: where in x ij , i represents each one of the measurements of j days after sowing,x j is the average for j days after sowing, and n is the total number of samples for the date. As for the interpretation this value, a low MAD value indicates that the data points are closely together; while a high MAD value indicates their dispersion.  Table 1 suggests that there is more dispersion (higher MAD values) in LAI measurements in the reproductive phase (for 67 and 82 days after sowing) with 1.44 and 1.51, respectively; than in the vegetative (for 47 days after sowing) with 1.12 and ripening phases (for 97 and 116 days after sowing) with 0.66 and 0.45, respectively.
In general, these variations in LAI measurements for the same collection day can be attributed to the nature of the experiment. As the experiments were carried out on a commercial farm, with particular soil and terrain characteristics. Basically, the terrain was not a flat and well-maintained crop surface, but an even surface where muddy puddles were formed, which led to variations in plant growth. Moreover, as in any commercial farm, the seed planting was performed via a broadcast seeder; therefore, the location of plants is not equally spaced, specially when compared to traditional experimental plots where row planting and drop spreaders/seeders are often used. Finally, it is known that the growth stages for IDIAP 52-05 work as a gradient with some plants changing phases more rapidly than others [33]. These growth variations are reflected in our LAI measurements.

Model for In-Field Spectral Estimation of LAI
Equations described in Section 2.7 were used to understand the resulting relationship between in-situ reflectance and True LAI in the vegetative and reproductive phases, with n = 52 points averaged into 4 points representing dates after sowing (47, 67, 82 and 97). Model #1, for NDVI, resulted in a significant linear model (p < 0.05) with a multiple Rsquared of 0.92 and an adjusted R-squared of 0.89, and a RMSE of 2.22. The resulting model can be expressed as LAI = 8.57 × NDVI. As for Model #2, for MTVI2, resulted in a statistically significant linear model (p < 0.05) with a multiple R-squared of 0.92 and an adjusted R-squared of 0.89, and a RMSE of 2.20. The resulting model can be expressed as LAI = 10.45 × MTVI2. These results show that the MTVI2 base model has a slightly higher predictive ability of LAI, with a slightly smaller error than the one using NDVI.  Table 2 presents the ratio of change between the correction factor against irradiance by satellite band measured with the field spectroradiometer. Table 3 shows correction factor regression with solar irradiance W m 2 per band was performed to determine the change in correction factor due to changes in solar irradiance.

Relative Correction of the Satellite Image
Correction reflectance applied according to Equation (11) earlier for the irradiance of 8 February 2018, suggests a correction of 0.61911 for the Blue Band. While corrections of 0.6121, 0.4058, and 1.355 must be applied to the Green, Red, and NIR Bands, respectively.
Each model was applied on a PlanetScope satellite image to obtain the LAI maps (see Figure 9). The image is from 8 February 2018 and the object of observation was the crop IDIAP FL 137-11.   We worked with the exact polygon measured in the field. For this purpose, the coordinates of the points measured in the field were used. The points of the edges and ends were taken in order to form a polygon. The coordinates of the polygon are 4 points as listed below:

Validating the Satellite Image Correction
Tables 4 and 5 shows the LAI prediction error for Model #1 and Model #2, respectively. For both models the estimated LAI with the PlanetScope bands against the true LAI (obtained by cutting leaves) are compared for both models. Estimation from Two LAI measurements was employed for each model. Model #1's average percentual error is 13.52%; for Model #2 the error amounted to 18.65%. Both cases required a detailed analysis of the relative correction of the satellite image. This analysis was made to understand further the model and how this error could be reduced.

Discussion and Conclusions
Visualizing spectral signatures serves to identify bands that can be used to discriminate between classes. Yang and Chen [47] mention that in the vegetative phase, the signatures are affected by: the change in leaf area, the growth of the leaves in number and size, and the soil discovered; the smaller the leaf area, the greater the interference of the soil discovered in the firm. The spectral response of rice cover crops grown within a few days after sowing showed a constant variation in the reflectance range along the signature. This spectral response does not reach the shape of a typical plant toppings signature. Spectral responses at ages 47, 67, 82, and 116 days after sowing show the typical distinctive features of vegetation: low reflectance in the visible light range, a noticeable peak in the green around 570 nm, and high near-infrared reflectivity [10].
With the plant's growth, an average increase in reflectivity is observed in the nearinfrared, but the variability in this band is high, suggesting difficulty in detectability. The variation in the near-infrared range is likely due to differences in the water content in the plant's leaves [10].
The reflectance was higher at 67-82 days after sowing (which comprises 2 phases, the final days of the reproductive phase and the beginning of the ripening phase) compared to the other sowing dates, which coincides with the results of [48], in tests performed on the Kongyu rice variety in China. The presence of the panicle also increases reflectance in the near-infrared range, which is consistent with the results of He et al. [30] for the japonica rice variety. At 116 days after sowing (last stages of cultivation), a decay of reflectance is observed, possibly due to decreased leaf water content and the loss of green leaf area.
Relative to the LAI, each true LAI value resulted from averaging three measures of LAI for each planting age. In the measurements of true LAI of the variety IDIAP 52-05 it is observed that the samples of an area of 25 × 25 cm 2 present high variability especially at the older age of the crop.
As for the treatment of the images, these are acquired with a correction that estimates the reflectance at the level of the Earth's surface [35]; which consists of the application of the standard atmospheric correction that is made with the algorithm 6S version 2.1, with continental aerosol content model and optical thickness data provided by MODIS satellite images [49]. The atmospheric correction methodology applied to PlanetScope products has some limitations such as needing overlap with areas measured with the MODIS satellite and the lack of correction in the presence of thin clouds or mist [49]. The additional correction factor brings the measurement of PlanetScope's corrected product at a ground level closer to the direct reflectivity measurement with the spectroradiometer. The correction factor is not static and varies between bands as it can be seen in Table 3. When comparing the Measured LAI for the verification points shown in Tables 4 and 5 to the average spectral signature of plants with the same days after sowing (with in-field measurements and those shown in Table 1) it seems plausible and very well in range with a plant with 77 days with 7.71 and 6.51, respectively, thus a verification of the procedure.
This work applied individual band correction factors to the satellite images. For all bands, a quadratic model was fitted. The correction factor tended to decrease to a minimum of irradiance defined for each band with values of 110 (blue), 185 (green), and 175 (red), and then increase. For the case of the NIR band, the behavior of the correction factor was precisely the opposite, showing an increase to a maximum irradiance of 160. Then a decrease in IR band. Table 2 shows the correction equations used for each band and the irradiance range applicable to the theme represented.
The validation test performed using two measurements of true LAI of the IDIAP FL 137-11 crop 77 days after sowing determined that model 1 presented an error rate of 13.52% comparing the estimated and true LAI values. Model 2 with MTVI2 obtained a predictive approximation with a 18.65% error. In general, literature shows that MTVI2 has a higher correlation rate with LAI for corn, soybeans, and wheat and less sensitivity to chlorophyll variations [13], such as rice in the ripening phases.
More importantly, the average magnitude of the errors for both Model #1 and Model #2 is around 15%, it can be attributed to many factors, but it can be accounted to the uneven distribution of LAI values. In fact, the maps shown on Figure 9 show that LAI estimates on the images form big homogeneous groups for Model #1 and small heterogeneous groups for Model #2. These notable distribution differences might be due to irregularities at the ground level, a byproduct of the experiment being performed on a commercial farm instead of an experimental plot. In addition, errors of estimated LAI around 20% are most likely due to other limitations of this study, such as the inevitable change in rice varieties, differences in the soil, and the seasonal difference between collection and validation dates. That are mostly noted by points having the same LAI, but changing significantly in the in-field spectral signatures in the NIR bands. In addition, different LAI values for similar NIR band values.
Despite limitations, and the fact that the resulting fitted models (and their results here presented), are only valid for PlanetScope satellite images, since they were adjusted to the bandwidths and sensitivity of PlanetScope Dove classic equipment. This work provides an important advance in precision agriculture, specifically in monitoring LAI for rice crops. The LAI is a physical parameter closely related to productivity and key to applying the various models employed, including the one related to MTVI2, allowing predictions consistent with the real LAI using PlanetScope satellite imagery.

Data Availability Statement:
In-field spectroscopical data collected is available upon request.