Crop Monitoring Strategy Based on Remote Sensing Data (Sentinel-2 and Planet), Study Case in a Rice Field after Applying Glycinebetaine

Abstract

World agriculture is facing a great challenge since it is necessary to find a sustainable way to increase food production. Current trends in advancing the agriculture sector are based on leveraging remote sensing technology and the use of biostimulants. However, the efficient implementation of both of these on a commercial scale for the purposes of productivity improvement remains a challenge. Thus, by proposing a crop monitoring strategy based on remote sensing data, this paper aims to verify and anticipate the impact of applying a Glycinebetaine biostimulant (GB) on the final yield. The study was carried out in a rice-producing area in Eastern Spain (Valencia) in 2021. GB was applied by drone 33 days after sowing (tillering phase). Phenology was monitored and crop production parameters were determined. Regarding satellite data, Sentinel-2 cloud-free images were obtained from sowing to harvest, using the bands at 10 m. Planet data were used to evaluate the results from Sentinel-2. The results show that GB applied 33 days after sowing improves both crop productive parameters and commercial yield (13.06% increase). The design of the proposed monitoring strategy was based on the dynamics and correlations between the visible (green and red) and NIR bands. The analysis showed differences when comparing the GB and control areas, and permitted the determination of the moment in which the effect of GB on yield (tillering and maturity) may be greater. In addition, an index was constructed to verify the crop monitoring strategy, its mathematical expression being: NCMI = (NIR − (red + green))/(NIR + red + green). Compared with the other VIs (NDVI, GNDVI and EVI2), the NCMI presents a greater sensitivity to changes in the green, red and NIR bands, a lower saturation phenomenon than NDVI and a better monitoring of rice phenology and management than GNDVI and EVI2. These results were evaluated with Planet images, obtaining similar results. In conclusion, in this study, we confirm the improvement in rice crop productivity by improving sustainable plant nutrition with the use of biostimulants and by increasing the components that define crop yield (productive tillers, spikelets and grains). Additionally, crop monitoring using remote sensing technology permits the anticipation and understanding of the productive behavior and the evolution of the phenological stages of the crop, in accordance with crop management.

1. Introduction

A strong agricultural sector is of paramount importance if the global population is to maintain its current economic and social progress. Currently, 8.9% of the world’s population suffer from hunger [1], a worrying number considering the increase in world population predicted by the UN in the coming decades [2]. This is a 21st century problem, with an added environmental component—the increase in world food production must be sustainable [3]. This calls into question the green revolution of the 20th century, to which a large amount of global pollution is attributed [4]. In this context, some authors have pointed out the need for a second green revolution, which is capable of strengthening world food security in a sustainable and environmentally friendly way [5,6]. Although scientific publications have reflected this need, exhibiting a greater presence of genetic and technological improvements [7], these advances are a long way from being manifested in the current data on the production of raw food materials [8]. A clear example is the cultivation of rice, which contributes 21% of the calories ingested in the world, and thus is the most important food crop [9]. According to FAO statistics, its world production has come to a standstill in many countries, highlighting how the major world powers (China, USA and the European Union) suffer from this phenomenon. China, despite being the world’s largest rice producer, is also the largest importer [10], with its imports increasing considerably in recent years [4]; the USA, meanwhile, has reduced its participation in the international market due to the protection policies of Asian countries [10]. For its part, the European Union (whose production is limited, representing 0.4% of world production), is the second largest importer in the world [10]. This trend serves as a warning that each country has to face up to the challenge outlined by the second green revolution, reaffirming the need to design and apply new cultivation methodologies for the purposes of yield optimization [3].

In agriculture, several factors lead to the instability of international markets, some of which transcend human control [11]. In this way, abiotic stresses can generate unexpected fluctuations in food production; the cultivation of rice, given its great demand for water and climatic requirements, does not escape this trouble [12]. Drought, salinity and extreme temperatures greatly condition the crop. Faced with this problem, the use of biostimulants is considered a valid alternative in the scenario of the second green revolution [13]. A widely accepted biostimulant is a compatible solute called glycinebetaine (GB). It is a fully N-methyl-substituted derivative of glycine, and it can be found naturally in many plants [14]. Its beneficial effect on abiotic stresses is widely recognized, and in some crops, it is even able to improve yield under non-stressful conditions [15]. Rice cultivation is not capable of synthesizing GB naturally, so exogenous applications considerably increase resistance to abiotic stresses [16,17]. The effect of GB on the crop is known at the physiological level, and the improvement in yield under stress conditions has been confirmed; however, there is limited information about crop behavior under natural conditions. Moreover, few studies have monitored the effect on crop growth from its application until harvest. In this context, it is difficult to find the determinant productive and physiological variables that affect the final yield, which highlights the need to introduce new variables into the agronomic evaluation.

The evaluation of the distribution of the productive parameters of the crop at within-field level has classically been carried out with discrete measurements in the field. Recently, with the consolidation of space missions for Earth observation (EO) in the solar spectral range (including visible (VIS) and near-infrared (NIR) spectra) at high spatial and temporal resolutions, it has been shown that this classic yield evaluation is not very effective, raising the need to incorporate remote sensing data to crop monitoring [18]. In other words, satellite data can be included to perform the agronomic evaluation of the effect of crop management on the final crop production, such as GB application. Novelli et al. (2019) [19] confirmed that yield forecast models can be improved by incorporating remote sensing data. With remote sensing, it is possible to identify the within-field level variability, which many crop models consider to be uniform [20]. However, the optimum integration of satellite data in order to monitor the whole set of crop agronomic variables that can, in turn, be used to improve crop management requires a huge dataset that is difficult to acquire (soil type, climate variables, management practices, phenology, plant genetics). Huang et al. (2018) [21] highlighted the need to simplify both the identification of the rice crop’s phenological stages and its management to improve the yield forecasting. In a wheat study, Skakun et al. (2019) [22] increased the determination coefficient between the spectral reflectances of satellite and marketable yield by incorporating a phenological adjustment using the accumulated temperature data from the crop. This method can potentially fit with temperature data when there is a significant difference between the studied regions. In Spain, there is great variability in the growth of the rice crop within a uniform climatic region. Franch et al. (2021) [23] showed that crop growth in fields cultivated on similar dates differs, which in turn affects the crop modelling of some phenological stages.

Since the launch of the first space mission to observe agriculture in 1972, there have been several studies looking for relationships between the spectral reflectance from satellite sensors and biophysical crop parameters [24]. In the VIS spectral region, more than 90% of solar radiation is absorbed by the leaf pigments in the blue and red bands, and more than 80% in the green band [25]. In fact, Thomas and Gausman [26] showed that these bands are correlated with the chlorophyll content, and that this increases in line with the reflectance decrease. Meanwhile, the NIR spectral region, where there is hardly any absorption of solar radiation, does not show any significant correlations with the pigments of the vegetation [24]. It is, however, tightly correlated with the crop biomass and leaf area (Leaf Area Index, LAI) [27], showing a greater reflectance in the NIR with higher biomass or LAI. Note that both parameters and the chlorophyll content are positively related [28]. Given these differences between the properties of the visible region and the NIR, a multitude of vegetation indices (VI) have been designed with the aim of integrating all the biophysical parameters of the crop in the same value, so the variations due to atmospheric effects are reduced [29]. Thus, one of the most widely used indices is the NDVI (normalized difference vegetation index) [30]. However, the saturation phenomenon of this index for high LAI values (>2) has been widely demonstrated [31,32,33], partly caused by the red band, which remains unchanged for higher chlorophyll and biomass content [33,34]. Meanwhile, the NIR band, with a high percentage of transmissivity, presents a lower saturation [35]. As a consequence, a multitude of corrections have been proposed by the NDVI to reduce the saturation effect [36], such as changing the red band for the green band, which presents a smaller saturation phenomenon [37], as well as the use of the red edge bands [38,39]; or using the difference vegetation index (DVI), based on the simple difference between the red and NIR bands, which is well-correlated with crop yields and responsive to high yield values [40,41]. Focusing on each band, the blue band presents problems in the interpretation of its values due to the marked influence of aerosols in this region [42,43]; on the other hand, at a theoretical level, the reflectance in this region is very similar to that occurring in the red region [37,44]. In addition, this band is less sensitive to chlorophyll content than the others [24], and while only chlorophyll absorbs radiation in the red region, absorption by carotenoids is also involved in the blue region [44], all of which justifies the difficulty of incorporating the blue band in the VI design. Of the two remaining bands in the visible range, red is the most commonly studied region; this band has a close relationship with the NIR reflectance [45,46]. On bare soil, some studies [47,48] have verified a close positive linear correlation between both bands. However, as the crop covers the entire surface, the reflectance in the red band tends to decrease, while the NIR’s reflectance increases. Most of these studies have been carried out on older satellites or with measurements taken in the field. Nowadays, with the existence of new satellites capable of providing data with a high degree of reliability and precision (for example, Sentinel-2 and Planet), and with considerably increased spatial, spectral and temporal resolutions, it is possible to assess the performance of these indices using these new platforms. Franch et al. (2021) [23] highlighted how, in rice cultivation, it is possible to use satellite data to monitor its final yields. Nevertheless, the spectral response of the crop could vary depending on its phenological stage and agronomic management [33,49]. Considering all of these parameters, the design of a crop monitoring strategy with satellite data throughout its entire cycle is still a challenge to overcome.

Through the use of remote sensing data, this paper aims to verify and to anticipate the productivity improvement in the JSendra rice variety by applying GB at the optimal time (tillering stage) [23]. This GB effect is intended to be monitored with Sentinel-2 and Planet data throughout its season, identifying its phenological stage and crop management. Two scientific questions are proposed as a means of achieving the main goal. Could a monitoring strategy be designed based on the evolution of the linear correlations between the reflectance of the red, green and NIR spectral bands for each date? Can a vegetation index (VI) maximize sensitivity to crop biophysical changes in the visible and NIR bands, integrating all of the changes into a single value? The objective of the VI study is to find the most suitable VI with which to identify the interactions between the visible and NIR bands in the rice crop after modifying either the agronomic management (GB application) or the plant phenology. The VIs used for crop monitoring have traditionally considered correlations with physiological parameters; the correlations of each VI with the spectral band response, crop management and phenology in rice have not yet been analyzed with the new higher resolution satellite missions (Sentinel-2 and Planet); thus, this study presents an in-depth analysis of the rice crop.

2. Materials and Methods

2.1. Study Area

The experiment was performed during 2021 in a traditional rice-producing area at sea level (Albufera wetland), in the coastal region of Valencia (Spain). The costal wetland of Albufera has an area of 211.2 km² and is bordered by the Turia and Jucar rivers located to the north and south, respectively. The Albufera is the second-largest lake in the Iberian Peninsula, with surface area of 23.2 km². The European Commission [50] restricts agricultural practices in the area to just rice crops since it is considered a special protected area in the Natura 2000 network. This area can be considered as a homogeneous rice planting area of approximately 10 × 20 km² in extension [23]. Figure 1 shows the location of the experiment in the region.

According to the Papadakis’ agroclimatic classification system [51], the climate in this area is subtropical Mediterranean with hot and dry summers. The soil has the following properties: sandy loam; pH: 7.9; organic matter: 3%; and electrical conductivity (EC): 3.3 dS·m⁻¹. The irrigation water comes from the Albufera lake, and has no rice crop restrictions regarding salinity (pH: 7.5; EC: 3.2 dS·m⁻¹), according to Ayer and Westcott (1994) [52]. Water management is carried out by flooding the fields (to a depth of 15 cm) during the whole season, with the exception of three periods, when the fields are traditionally dried for terrestrial treatments over an entire irrigation sector.

2.2. Experimental Design

JSendra, a common japonica-type Spanish rice cultivar, was used in the experiment. Rice was sown on 24 May (0 DAS, days after sowing) in Silla (Valencia, Spain), using a sowing rate of 180 kg·ha⁻¹, and the harvest took place on 12 October (141 DAS). Nitrogen, phosphorus and potassium were applied as described by Osca (2013) [53], and the rates were as follows: 170 kg N·ha⁻¹; 60 kg P₂O₅·ha⁻¹; 150 kg K₂O·ha⁻¹.

The experiment consisted of four replications of 5000 m², each one conducted following a completely randomized design (Figure 1). One treatment, with two concentrations: 0 and 5.0 L·ha⁻¹ (1450 mM) of GB biostimulant was studied. Note that hereafter the statistics representing the control and the GB correspond to the average of the four replication areas in each case. GB was applied on 33 DAS (tillering stage), according to the rice crop model obtained in Franch et al. (2021) [23]. The GB was applied by drone to avoid any additional machinery overpass: AGR model A6; tank volume: 6 L; flow rate: 0.57 L·min⁻¹ by nozzle; pressure: 6 bar; forward speed: 6 m·s⁻¹; spray height: 1.5 m; swath width: 2 m; and 4 nozzles (Agroplast Ppij, model 6MS 01C anti-drift orange with ceramic tip, fan-shaped, 110° [54,55]).

Phenological stages were classified according to the BBCH scale [56]. Figure 2 shows the phenological cycle of the crop in the experimental plot and the periods when the fields are dried. Climatic data were registered by a datalogger [57]. Figure 3 shows the most important meteorological crop growth parameters that were retrieved: air temperature (T) and relative humidity (RH) (no rainfall occurred during the drying periods).

2.3. Determination of Production Parameters

Plant height, panicle length, panicles per m², filled grain per panicle, weight of 1000 grains, grain length and width were obtained by means of field measurements at harvest time. These measurements were acquired in an area of 0.25 m², randomly selected avoiding the edge effect (10 m inside) in 4 repetitions for each replication. In addition, we measured the days elapsed from sowing until the emergence of 50% of the panicles and the marketable yield (kg·ha⁻¹) for each replication.

2.4. Satellite Data

Satellite images obtained by the Multi-Spectral Instrument (MSI) on board the Sentinel-2A/B constellation [58] of tile T30SYJ were used. Sentinel-2 (S2) captures images of the Earth’s surface in 13 spectral bands; from these, the bands with the highest spatial resolution (10 m) were considered in this study: blue (central wavelength ≈ 490 nm), green (≈560 nm), red (≈665 nm) and near-infrared (NIR) (≈842 nm). Cloud-free images from sowing to harvest in 2021 and an image from 20 April 2021 were selected prior to sowing. The downloaded product was level 2A, which provides surface reflectance [59].

Images were also obtained from Planet satellites belonging to the private company Planet Labs [60]. From the Planet image catalogue, only cloud-free PS2 images were used in this study from the date when the crop predominates over the soil effect until harvest. Data were acquired from the NASA program [61]. Planet constellation presents 4 spectral bands with a spatial resolution of 3 m: blue (central wavelength ≈ 485 nm), green (≈545 nm), red (≈630 nm) and near-infrared (NIR) (≈820 nm).

Table 1. Dates studied for each satellite, with DAS standing for days after sowing.

Sentinel-2		Planet
DATE	DAS	DATE	DAS
20 April 2021	−34	8 July 2021	45
30 May 2021	6	13 July 2021	50
9 June 2021	16	29 July 2021	66
14 June 2021	21	7 August 2021	75
24 June 2021	31	25 August 2021	93
4 July 2021	41	26 September 2021	125
14 July 2021	51
19 July 2021	56
29 July 2021	66
28 August 2021	96
2 September 2021	101
12 September 2021	111
17 September 2021	116
7 October 2021	136
12 October 2021	141

shows the dates studied for each satellite.

2.5. Methods

Figure 4 shows the workflow proposed in this article to monitor the crop with remote sensing data. Each analysis is related to the identification of a variable or phenomenon needed for crop monitoring.

2.5.1. Sentinel-2 Data Analysis

In order to analyze the satellite data, the dynamics between the reflectance in the visible and NIR bands were first checked. Thus, with the aim of determining the moment when the rice crop covers the entire area, the red–NIR dynamics were analyzed for each available date from sowing to harvest.

The next step consists of obtaining the linear correlations between both bands for all the available dates, adding the linear correlations between the green and the NIR bands to the analysis, as well as the green and the red bands, given their importance in crop monitoring [37]. Additionally, we added a satellite image taken one month before sowing to the analysis, a time when the soil was bare and completely dry, to contrast with the crop canopy growth with the bare soil when both flooded and dry.

Three widely used vegetation indices (VIs) were compared, namely the NDVI [30], GNDVI (green normalized difference vegetation index) [37]) and EVI2 (enhanced vegetation index [36]), and a new index was designed. The equations of NDVI, GNDVI and EVI2 are:

$$NDVI=\frac{NIR-Red}{NIR+Red}$$

(1)

$$GNDVI=\frac{NIR-Green}{NIR+Green}$$

(2)

$$EVI2=2.5\frac{NIR-Red}{NIR+2.4·Red+1}$$

(3)

where all the variables represent spectral surface reflectance.

Considering the spectral properties of the vegetation in each wavelength, we analyzed the linear correlation between each index and the spectral bands for the dates when the crop covers the entire surface [24,62].

The NDVI evolution was used as an indicator of the phenological stage of the crop, as suggested by Mosleh et al. (2015) [63].

Finally, a comparison of the averages of each VI across the control treatment and the GB treatment area was made on the selected dates.

2.5.2. Planet Data Analysis

Firstly, we analyzed an image from Planet satellites on the same date as Sentinel-2 (29 July (66 DAS)). Planet satellites acquired the data at 11 h 17 min, while S2 acquired them at 10 h 46 min. The short time elapsed in data acquisition between these two satellites makes their comparison possible. Planet data were aggregated to the 10 m resolution of S2, obtaining the r² coefficient of determination between each Planet band and the corresponding S2 band. In addition, the r² between the VIs calculated with Planet and the values obtained with S2 were calculated. Given the robustness of the S2 measurements and the known calibration problems for the Planet constellation, these correlations are useful for determining the best VI to minimize the differences between both satellites [64].

Secondly, the available Planet dates at 3 m spatial resolution, during the period when the crop covered the soil surface, were considered so as to compare the control and treated areas, and these were intercompared with the S2 results.

2.5.3. Software and Statistics

Statistical analyses were performed with StatGraphics Centurion 18 software [65]. The processing of the satellite images was carried out with the QGIS 3.10.14 software [66]. A comparison of means between GB treatment and control on productive parameters, spectral band reflectance and VIs were analyzed using the least significance differences statistical test (p < 0.05).

3. Results

3.1. Productive Parameters Analysis

Table 2. Influence of glycinebetaine (GB) treatments on productive parameters of rice plants.

Parameters	GB Plants	Control Plants	p Value
Plant height (cm)	80.1	79.2	ns
Panicle length (cm)	16.7	15.6	ns
Days to 50% panicle emergence	75	74	ns
Panicle per m2	321 a	305 b	**
Filled grain per panicle	96 a	92 b	**
Weight of 1000 grains (g)	35.8 a	34.5 b	**
Grain length (mm)	8.10	8.05	ns
Grain width (mm)	4.12	4.04	ns
Grain yield (kg·ha−1)	10767 a	9523 b	**

** and ns: significant at p < 0.01 and non-significant. Different letters for each parameter indicate statistically significant differences using LSD test (p < 0.05).

shows the effect of the treatments (GB and control) on the productive parameters of the crop. The exogenous application of GB does not show a statistically significant influence on the morphological parameters of the plants and grains, or on the moment of emergence of 50% of the panicles. Nevertheless, in the three components of the yield (panicle per m², filled grain per panicle and weight of 1000 grains), the differences between treatments are statistically significant (p < 0.05); each value is higher in GB plants. These differences lead to a higher marketable yield in plants treated with GB (p < 0.05) (increase of 13.06%).

3.2. Sentinel-2 Data Analysis

3.2.1. Dynamics and Correlations between Visible and NIR Regions

Figure 5 shows the dynamics between the reflectance in the red and NIR bands for the available dates from sowing to harvest. From 6 DAS to 31 DAS, the reflectance in the red band does not decrease. At 21 DAS, the fields were dried for the purposes of chemically controlling the weeds competing with the crop, generating a considerable increase in band reflectance. At 31 DAS, the fields were flooded again, reducing the reflectance, demonstrating that water management exerts an influence on it. As the proportion of vegetation increases, the NIR bands also increase their reflectance. However, as with the red band, water management has a notable influence: at 21 DAS, an increase in reflectance is generated that is greater than that of the crop at 31 DAS, with statistical differences at p < 0.05 (

Table A1. Table of means of reflectances in the red and NIR bands over the whole crop cycle, combining control and GB. Different letters for each band indicate statistically significant differences between dates using the LSD test (p < 0.05).

DAS	Red Reflectance	NIR Reflectance
6	0.0597 a	0.0970 a
16	0.0752 b	0.1204 b
21	0.1129 c	0.2068 c
31	0.0658 d	0.1434 d
41	0.0330 e	0.1959 e
51	0.0287 f	0.3405 f
56	0.0164 g	0.2809 g
66	0.0157 g	0.3414 f
96	0.0198 h	0.3702 h
101	0.0222 i	0.4080 i
111	0.0402 j	0.4015 j
116	0.0426 k	0.4047 ij
136	0.0710 l	0.3048 k
141	0.0766 b	0.2969 l

).

Although a decrease in red reflectance is observed at 31 DAS, it is from 41 DAS onwards that it decreases significantly (p < 0.05, Table A1) with respect to that existing at 6 DAS. These differences are maintained until plant senescence (116 DAS, Table A1).

At 51 DAS, the fields are dried according to the agronomic management of the crop (to enhance tillering). At this moment, the reflectance in the NIR increases because of the elimination of the water sheet, which is coherent with the Kimura et al. (2004) [47] findings. However, in this second drying, the reflectance in the red band does not increase, but rather slightly decreases (p < 0.05, Table A1). Once the fields are flooded again (at 56 DAS), the reflectance values in the NIR decrease due to the background effect of the water sheet. The NDVI peak for the available dates occurs at 66 DAS (Figure A1), which coincides with the panicle initiation stage.

From 96 DAS onwards, the crop is in the ripening phase, and at 136 DAS, it starts senescence. During these dates, the reflectance in the red increases, reaching the maximum value at harvest (141 DAS), showing statistically significant differences with respect to the rest of the dates (Table A1); in the NIR meanwhile, there are no statistically significant differences between dates (Table A1), and it decreases on the dates close to harvest (at 136 and 141 DAS).

Figure 5 also shows a great variability in the relationships between the red spectrum and the NIR. In fact, Figure 6 shows the evolution of these correlations between different bands (NIR vs. red; NIR vs. green; red vs. green) on every date under consideration. The dry bare soil (−34 DAS = 34 days before sowing) and flooded bare soil (6 DAS) show a positive correlation between the bands. As the growth of the crops begins (from 16 DAS onwards), the correlation coefficient between red and NIR decreases and is not statistically significant up to 41 DAS. From 41 DAS to harvest, the correlation coefficient between these bands is negative and it is again statistically significant (p < 0.01), reaching the minimum values at 41 and 51 DAS. From these dates, a saturation of the surface reflectances and, consequently, of the NDVI, is verified (Figure A1). At 96 and 101 DAS, new minimum values of the correlation coefficient between NIR and red are reached. The correlation between green and NIR presents a greater variability than that obtained in the previous correlation.

On the other hand, the correlations between green and red fluctuate in a range of r of between 0.4 and 0.93. These values are statistically significant (p < 0.01). The maximum r values for this correlation are reached at 6, 41, 51 and 136 DAS, and they coincide with the minimum of the correlations between the NIR and the two visible bands, except for 6 DAS. During the mature phenological stage, it may be observed that the correlation between the two visible bands increases, while the correlations between the NIR and the visible bands decrease. In addition, it can be seen that the correlation coefficient between NIR and red reaches the highest absolute values compared to the correlation between NIR and green.

According to the results of Figure 5 and Figure 6, we consider that, at 41 DAS, the crop predominates over the soil (the reflectance in the red band decreases significantly and the correlations between NIR and the visible bands are negative and significant).

3.2.2. Effect of GB on Spectral Band Reflectances

Figure 5 and Figure 6, show that on the date chosen for the application of the biostimulant (33 DAS), the crop still does not cover most of the soil surface, so the soil effect predominates over the plant (initial phase of tillering). Before treatment (31 DAS), no statistically significant differences were observed between the control and the GB treatment of the reflectance in the visible (red and green) and NIR bands (data not shown). On the first date available after the treatment (41 DAS), the analysis of the red–NIR dynamics (Figure 5 and Figure 6), outlined in the previous section, shows how the effect of the crop already predominates on the soil.

Figure 7 shows the evolution of the reflectance in each band for the treated crop surface (GB) and the untreated surface (control). We divided the analysis into dates before (before NDVI peak) and after (after NDVI peak) flowering, given the morphological and growth changes that occur [63]. Statistically significant differences (p < 0.05) were detected between the plants treated with GB and the untreated on every date in the green and NIR bands. However, the red band reflectance does not show statistically significant differences either during the panicle initiation (66 DAS) or the beginning of ripening (96 DAS), with differences appearing in the vegetative phase and in maturation and senescence. The statistically significant differences between the treated and untreated plants during the maturation and senescence phase in the three bands remain for up to 108 days after the GB treatment (harvest). During the mature stage of the grain and the senescence of the plant, the evolution of the values in the three bands is uneven: in the green band, the values stabilize; in the red band, there is a notable increase in reflectance; and in the NIR band, there is a significant decrease after a remarkable plateau.

3.2.3. Construction of the Vegetation Index

The difficulty in correlating the VI with the biophysical and productive properties of a crop is well-known, and for this reason, various authors propose the use of alternative indices that incorporate three spectral bands with the aim of increasing the sensitivity to changes in the crop canopy [67]. These studies have mostly been carried out with field measurements (spectroradiometer), so that the spectral bands of red edge, SWIR (short-wave infrared) and blue bands are used at very high spectral and spatial resolutions, without the need for atmospheric correction.

Currently, satellite missions at an adequate spatial resolution for existing rice fields in the study area (10 m) only have three bands in the visible and one in the NIR range. On the other hand, Franch et al. (2021) [23] have shown how the incorporation of the red edge and SWIR bands in the modelling of the crop in the study area barely improves the results. It has been found by other authors that the greater the biomass, the more reflectance there is in the NIR, while in the visible bands, the reflectance decreases after such an increase [28]. Therefore, the result of the difference between the reflectance in the NIR and the reflectance in the visible range can be indicative of the state of the vegetation cover, this principle being the same one used in the formulation of the NDVI and GNDVI [28,37].

By taking advantage of the influence of vegetation on the NIR spectral range, we propose the use of an index that integrates the NIR and visible bands in order to maintain a high spatial resolution (generic equation). However, due to the greater atmospheric influence in the blue band [36,43,68], it was not considered in the generic equation.

The generic equation for the construction of the index:

$$Index=\frac{NIR-\sum VIS}{NIR+\sum VIS}$$

(4)

The equation of the designed index, called “Normalized Crop Management Index” (NCMI):

$$NCMI=\frac{NIR-\left(Green+Red\right)}{NIR+\left(Green+Red\right)}$$

(5)

where all of the variables represent spectral surface reflectance.

Table 3. Average of r between each band and VI by performing a linear regression for the available dates (>41 DAS) on Sentinel-2 (the last row is the average of the r of every band in terms of absolute value).

	NDVI	GNDVI	EVI2	NCMI
NIR	0.82	0.80	0.94	0.83
Red	−0.96	−0.81	−0.74	−0.90
Green	−0.70	−0.90	−0.59	−0.83
Average	0.82	0.83	0.76	0.85

shows the mean value of the correlation coefficient between each spectral band and VI for the dates analyzed in this study (>41 DAS). From the results obtained, it may be seen that the correlation coefficient of the reflectance values in the NIR is higher with EVI2; that of the reflectance of red is greater with NDVI; and that of the reflectance of green with GNDVI. Individually, the NCMI does not improve any previously noted correlations; nevertheless, in the global assessment of the NIR, the red and the green bands, the mean correlation coefficient (r = 0.85) in this VI is the highest compared to the rest of the VI.

Figure 8 shows the evolution of the differences between the average value of each VI in the crop area treated with GB and that used as control from 41 DAS until harvest. There are fewer differences as the crop approaches the NDVI peak (maximum greenness [63]); however, once the crop is in full maturity, the differences increase again, so that there is a relative maximum of the difference at 136 DAS. The differences between treatments (GB and control) are statistically significant for every VI on each date, except for NDVI at 66 and 96 DAS (Figure A1). NDVI is the VI with the smallest differences between treatments, followed by GNDVI.

NCMI and EVI2 are the indices that show the greatest differences, with NCMI showing greater differences on every date, except at 96 DAS. Finally, NCMI and EVI2 are the VIs that most closely resemble the achieved percentage yield increase (13.06%), with NCMI performing better at later dates, especially at 136 DAS (fully formed grain).

Figure 9 shows the NCMI values for each date analyzed, confirming the differences between treatments on the aforementioned dates, the phenology and the spatial distribution across each region.

3.3. Planet Data Analysis

Table 4. Coefficient of determination (r²) between the data of Sentinel-2 and Planet at 10 m spatial resolution for each band and VI at 66 DAS (p < 0.01 for each band and VI except for the blue).

	r2
Blue	0.13
Green	0.56
Red	0.63
NIR	0.79
NCMI	0.77
NDVI	0.72
GNDVI	0.72
EVI2	0.80

shows the coefficient of determination (r²) obtained by correlating each S2 band and VI with its corresponding Planet value aggregated to 10 m for the only coincident date of satellite data acquisition (29 July 2021, corresponding to 66 DAS). The results show that the Planet and S2 correlation depends on the spectral band analyzed. The blue band does not present any correlation, while the rest of the bands (green, red and NIR) present a significant correlation (p < 0.01), the value of r² rising as the wavelength increases. Regarding VIs, EVI2 shows the best correlation (r² = 0.80), with NDVI and GNDVI performing the worst. Aside from being the most sensitive VI to changes in the visible, even when using Planet data (data not shown), NCMI presents an r² close to that of EVI2 and to the NIR band.

Figure 10 is analogous to Figure 8, and compares the two treatments based on all available Planet images starting from the date when the crop signal predominates over the soil (>45 DAS). The evolution of all the indices is very similar to that obtained with S2 data. The NDVI peak for the available Planet dates occurs at 75 DAS (data not shown). NCMI is again the index that shows the greatest differences between treatments (GB and control) and the index that determines the percentage difference that is closest to the percentage increase in yield (13.06%). The results obtained from Planet clearly show two groups of VIs: NCMI and EVI2 on the one hand and NDVI and GNDVI on the other. The differences between treatments are statistically significant for every date and VI (p < 0.05, data not shown).

Figure 11 shows the NCMI values at 3 m spatial resolution for each date analyzed. Similar to Sentinel-2 data, the differences between treatments and the phenological evolution of the crop are verified. The higher spatial resolution of Planet increases the analyzed area, reducing the edge effect. The availability of Planet images permits the identification of the evolution of the crop on dates not available in Sentinel-2, showing that the highest NCMI values are obtained at 75 DAS; however, the results show a constant bias compared to Sentinel-2. Figure A1 (Appendix A) shows these differences, and the behavior of each VI can be visually observed. The differences in values between Planet and S2 occur for every VI studied, and the Planet values are always lower. This effect might be caused by calibration issues of the Planet visible bands, as highlighted by previous authors [18].

4. Discussion

In Mediterranean climate conditions, rice is cultivated in flooded fields due to its water demands, and the crop cycle is concentrated during the hottest months of the year due to the temperature and photoperiod requirements. Consequently, high temperatures and a low relative humidity can condition the potential yield of the crop. Although the mean values are adapted to the physiological optimum of the crop [69], the extreme values may differ greatly from such mean values in critical phases, with high temperatures (T > 33 °C) coinciding with a low relative humidity (RH < 40%) in every stage (Figure 2). Therefore, under these conditions, the crops can be subjected to intermittent and random stresses during their cycle [70,71].

GB is identified in the scientific literature as an effective compound for combating abiotic stresses [14]. Generally, these stresses are induced by humans under very controlled conditions, so the effect of GB in the field under natural stress, that is, not forced, is not usually considered in studies. In this study, the positive effect of the exogenous application of GB on the rice yield components and its marketable yield is verified. This effect coincides with that found in previous studies by other authors [17,72]. The improvement in the productive potential of the crop is also observed in transgenic rice plants capable of accumulating the solute [73].

JSendra is a rice variety that has been improved with the objective of maximizing the harvest index (filled grain dry matter divided by total plant dry matter) and remote sensing technology allowed us to highlight that its yield depends strongly on the number of productive tillers per plant [74,75]. When studying this variety, Franch et al. (2021) [23] verified the importance of the tillering phase for achieving high yields, obtaining the best correlation between reflectance in the NIR region and final yields during tillering. NIR is closely related to the amount of biomass and LAI [28,76], and a greater reflectance in this band may be indicative of more vigorous tillering. Therefore, early tillering with GB may have stimulated a better growth of productive tillers per plant, as presented in Table 1. The number of productive tillers is the component of the final yield that is defined in the vegetative phase. Despite the excellent meteorological conditions, with optimal temperature and RH conditions from the application of GB until the first satellite image (41 DAS), the reflectance in the visible and NIR ranges presented significant inter-treatment differences. This phenomenon may indicate that GB is also effective under non-stressful conditions; diverse authors performed laboratory experiments and verified this effect in the cultivation of both rice and in that of other crops, such as tomato, a crop which does not synthesize GB either [15,72]. In this sense, Chen and Murata (2011) [73] highlighted the need to carry out a field validation of the possible effect of GB under natural conditions (optimal and non-optimal), pointing out how, in both scenarios, this biostimulant could increase world production in the future.

Natural stress conditions occurred during the crop season. The first period of high temperatures (38 °C) and low RH (16%) occurred at 41 DAS, accompanied by a similar period at 43 DAS, which can induce stress on the plantation. Later, at 49 DAS, the fields were left to dry up as a result of the management carried out in the area. This coincides with another period of very adverse temperatures and relative humidity (Figure 2). In this study, we show that remote sensing technology can assist the agronomical interpretation through different crop stages across the rice season. When S2 and Planet acquire images at 51 and 52 DAS, respectively, the fields have remained dry for 2 and 3 days. At this time, the green band shows the maximum difference between treatments of the entire crop season, while the red band also maintains a considerable difference. As the NDVI peak approaches, the differences between treatments for bands and VIs are reduced. This trend could be attributed to the VI saturation phenomenon [33]. In the remote sensing discipline, this phenomenon has been considered a limiting one for data applicability. However, saturation is a natural phenomenon and can be explained from an agronomic point of view. In all vegetation cover, the photosynthetically active radiation (PAR) or absorbed visible light reaches a saturation point in which an increase in LAI or chlorophyll content does not affect a greater absorption [77,78]; according to other authors [79,80], the plant could reach the highest photosynthetic rate of the entire cycle. Thus, it is not surprising that there are fewer differences between treatments during this period This phase prior to the NDVI peak coincides with stem elongation and panicle initiation, a moment that various authors identify as the most optimal for an anti-stress treatment, given the sensitivity of rice to spikelet sterility [81]. However, our remote sensing monitoring shows how there are fewer differences between treatments during this period. In other words, the first component of yield (number of productive tillers) has already been defined; thus, an anti-stress treatment at this moment will only improve the fertility of the developed panicles. Kim et al. (2021) [82] show how high temperatures during the vegetative phase affect a lower number of panicles per unit area. This decrease is offset by a greater number of spikelets per panicle. However, Ono et al. (2013) [83] indicate that the final weight of the grain depends on the transpiration and the carbohydrates synthesized in the flag leaf; therefore, a lower number of spikelets in the panicle can improve the weight of the final grain, given that there is less competition [84]. Therefore, a greater number of panicles allows for a more efficient use of the energy intercepted by the plant and translates this into a higher yield, as confirmed in our experiment. Our results are consistent with those reported by these authors, verifying a greater number of filled grains in the plants that were treated with GB (Table 1).

After the second component of the yield (number of grains per panicle) has been defined, the crop begins the ripening phase (grain weight will be defined in this phase). Kim et al. (2011) [85] indicate that the weight of the grain depends on the duration of this phase: the longer the phase, the greater the weight gain. Despite reducing the inter-treatment variations in the phase of stem elongation, when maturity is already developing (from the 96 DAS image), the differences increase again for all bands and VIs. Note that the differences in the red band are consistent, indicating greater light absorption, which may be related to greater photosynthetic activity in the treated area [24,76,79]. In this way, satellite monitoring verifies that senescence occurs later in the treated area, maintaining its greenness for a longer period. This result is extremely important because it indicates that the effect of the GB is maintained until the harvest date. In addition, the fact that the differences in the dates near the NDVI peak were either non-existent or reduced shows that the greater number of differences after the peak is not a direct consequence of the effect of the GB on the tillering; in other words, the plants can stay green for a longer time. This phenomenon is supported by recent studies that indicate how GB can induce the expression of specific genes that increase tolerance to oxidative stress, thus improving photosynthetic activity [73,86].

To summarize, it can be confirmed that the green and red bands behave similarly in relation to the phenological response of the crop. However, the green band presents a lower saturation phenomenon, and, in the maturation phase, its temporal evolution changes considerably. While the reflectance in the green band increases and stabilizes quickly (111 DAS), the reflectance in the red band takes longer to increase and does so significantly as the crop approaches senescence: at the dates closest to harvest it presents a greater reflectance than the green region (degradation of chlorophyll and prevalence of carotenoids [26]).

Additionally, the reflectance peak in the NIR appears after the NDVI peak. This NDVI peak coincides with the minimum reflectance in the visible range, since the crop loses its greenness when the value of NDVI decreases, according to Mosleh et al. (2015) [63]. Our results indicate that this may be the consequence of an increase in reflectance in the visible; an increase that Thomas and Gausman (1977) [26] linked to the beginning of the degradation of chlorophylls. Therefore, NIR would not be a good indicator of changes in crop physiology during early maturation. Furthermore, the reflectance in the NIR hardly varies on 101, 111 and 116 DAS, dates on which there are highly significant changes in the visible region. It should be noted that the reflectance in the NIR depends on anatomical factors in the leaf, the most important being the multiple scattering of solar radiation in the intercellular spaces of the mesophyll [44]. This multiple scattering is responsible for the reflectance in this region, increasing as these spaces grow larger, and a higher LAI is found [44,87]. With maturity, there are many phenomena that modify the internal anatomy and physiology of the leaf, some of which could mask others; for example, as the leaves dehydrate, reflectance can increase, while the degradation of chlorophylls can cause this to decrease [87]. Despite this complex interpretation of the evolution of the NIR after the NDVI peak, this band exhibits great inter-treatment variation. This relationship between NIR and the final rice yield is coherent with that found by Franch et al. (2021) [23].

For optimal monitoring, this article shows how the evolution of reflectance in the visible is a better indication of the phenology of the crop. In this way, there is great difficulty in correctly interpreting the NDVI, since there is a high saturation of the index under conditions of high vegetation cover [33] and a low correlation with the Green band, as shown in Table 2. Thus, it is necessary to find alternative VIs [32,33,34]. In this study, we test a new, simple normalized VI, which integrates the properties of each visible band. The proposed index modifies the NDVI and GNDVI, discounting a higher proportion of radiation visible to NIR. From a mathematical point of view, by incorporating the third term, we notice that the total size (NIR plus red plus green) increases and the remaining amount NIR–(red plus green) decreases. Thus, the relative size decreases. This is also possible since the red and green bands belong to the same order of magnitude (10⁻²), which is strictly smaller than the order of NIR (10⁻¹). Our results demonstrate that NCMI enhances the differences between the two rice treatments (GB and control plants) over the entire crop cycle. Meanwhile, the NDVI does not highlight inter-treatment variations studied during the phases of panicle initiation and the beginning of ripening (decisive phases in yield monitoring). On the other hand, the GNDVI, despite the similarity of its performance to that of the NCMI, does not consider the red band and, therefore, does not offer an accurate interpretation in the ripening and senescence phase of rice cultivation (Figure A1). Finally, the EVI2 differs regarding the evolution between the VI and the maximum light absorption of the canopy, because its performance is very similar to that of reflectance in the NIR (Figure A2).

NCMI also performs better when correlating S2 data with Planet, a result that could improve the fusion of timeseries from both satellites [64]. Finally, with Planet, the results obtained in the monitoring are consistent with those found with S2, which provides validity and consistency to the proposed monitoring and to the newly designed index.

The presence of clouds in the main rice-producing areas of the world makes it difficult to design a monitoring strategy for the crop [4,63]. Fortunately, Eastern Spain is characterized by a lack of cloud cover; however, clouds are still an unavoidable phenomenon. According to the modeling proposed by Franch et al. (2021) [23], the unavailability of dates in the reproductive phase (from 66 to 96 DAS) should not affect the optimal monitoring of the GB effect. In this phase, this study found that remote sensing data provide little information on the final yield in fields without agronomic problems. Despite the lack of images during this phase, our article shows that, both before and just after flowering, the differences between treatments are minimal. In this phenological phase, moreover, other authors show how the crop stops growing, maintaining its photosynthetic rate and its greenness [63,80]; there is no change in the leaves, so the spectral response of the canopy also hardly changes (in our paper, the NIR reflectance shows a very slight increase between 66 and 96 DAS compared to previous dates). Planet data in Figure 7 verify our theory. On 75 DAS, when the panicles emerge, the differences between treatments are minimal and only increase with the arrival of maturity.

5. Conclusions

The results obtained in the present paper demonstrate how useful remote sensing data are for monitoring the productivity of the rice crop. Glycinebetaine is widely recognized as a potent biostimulant for crop yield improvement; however, its agronomic influence on the productive behavior of the rice crop and the monitoring of its effect under natural field conditions remains unknown. The application of GB at the beginning of tillering (the critical phase, which is determinant in the final yield) has led to a better performance of yield components, thus achieving a significant increase in the final yield. The results also verify the need to analyze the dynamics of reflectance in all spectral regions, maximizing the sensitivity to changes in any band with a new index. This index should be tested in more areas and for longer periods in order to offer clear evidence that it contributes to the better monitoring of the rice crop; therefore, a preliminary study is presented in this paper. The analysis of band dynamics and correlations permits the identification of crop phenology and management, quantifying the effect of GB on productivity with remote sensing data. Knowing when the GB has a significant effect on the crop allows it to be in the best condition in the event of a specific stressful event, such as high temperatures, low relative humidity, drought, pests and diseases, etc., so that production inputs can be optimized.

The results reaffirm the modelling of Franch et al. (2021) [23], showing how the end of tillering and the transition from ripening to senescence are the key moments for yield monitoring. Thus, this study adds to the knowledge required for developing predictive models for rice crops.

We consider that our article can help to face up to the challenge represented by the second green revolution, as we introduce a new methodology for a more accurate control of production inputs. This new methodology effectively links crop management with Earth observation, opening the door to precise, sustainable agriculture.