Mathematical Integration of Remotely-Sensed Information into a Crop Modelling Process for Mapping Crop Productivity

Abstract

Remote sensing is a useful technique to determine spatial variations in crop growth while crop modelling can reproduce temporal changes in crop growth. In this study, we formulated a hybrid system of remote sensing and crop modelling based on a random-effect model and the empirical Bayesian approach for parameter estimation. Moreover, the relationship between the reflectance and the leaf area index was incorporated into the statistical model. Plant growth and ground-based canopy reflectance data of paddy rice were measured at three study sites in South Korea. Spatiotemporal vegetation indices were processed using remotely-sensed data from the RapidEye satellite and the Communication Ocean and Meteorological Satellite (COMS). Solar insulation data were obtained from the Meteorological Imager (MI) sensor of the COMS. Reanalysis of air temperature data was collected from the Korea Local Analysis and Prediction System (KLAPS). We report on a statistical hybrid approach of crop modelling and remote sensing and a method to project spatiotemporal crop growth information. Our study results show that the crop growth values predicted using the hybrid scheme were in statistically acceptable agreement with the corresponding measurements. Simulated yields were not significantly different from the measured yields at p = 0.883 in calibration and p = 0.839 in validation, according to two-sample t tests. In a geospatial simulation of yield, no significant difference was found between the simulated and observed mean value at p = 0.392 based on a two-sample t test as well. The fabricated approach allows us to monitor crop growth information and estimate crop-modelling processes using remote sensing data from various platforms and optical sensors with different ground resolutions.

1. Introduction

Satellite-based remote sensing is a useful technique to acquire spatiotemporal data consisting of a large number of pixels, but a relatively small number of temporal data points, from an agricultural field. Remote sensing can provide an inexpensive and non-destructive method from various platforms to collect a vast amount of information from an agricultural field [1]. However, there are two main difficulties related to the practical use of remote sensing data from satellite-aboard sensors in monitoring crop growth. Acquisition of remotely sensed data with a high spatial resolution that is feasible for use is limited, due to dependence on revisiting times, operational schedules, and clear sky conditions in the case of the optical satellite-based remote sensing that is more frequently applied in agriculture. Leaf area index (LAI) is commonly used to describe the canopy growth of a crop, but it cannot be measured directly from remote-sensing data. These contain the reflectances of light spectra across all waveband regions to measure the index of “greenness” of the plant canopy [1]. Meanwhile, appropriate crop growth models can provide temporal descriptions of crop conditions during the growing season, although the accuracy of this information is based on the quality of the input data, and on model design [2]. By effectively combining the advantages of remote sensing and crop modelling, the strengths of each approach may make up for inherent weaknesses in individual strategies [3].

There have been earlier efforts to combine the techniques of crop modelling and remote sensing. Arkin, et al. [4] proposed the concept of a hybrid model able to use Landsat data. This concept was included in grain sorghum growth simulation models, the sorghum growth model with feedback capacity (SORGF) [5] and the grain sorghum crop growth model (SORKAM) [6]. These models allow a parameter affecting leaf expansion rate to be adjusted to improve the agreement between simulated and measured LAI. Barnes, et al. [7] modified a cereal crop simulation model, CERES-Wheat [8] to allow the model to accept observed LAI values and to adjust related parameters in the model as a function of LAI. Similar approaches have been made using the WOrld FOod STudies (WOFORST) model [9] and the Simple and Universal CROp growth Simulator (SUCROS) model [10] to improve the overall model performance [11,12,13]. WOFORST is a simulation model for the quantitative analysis of the growth and production of annual field crops. The SUCROSE model simulates both potential and water-limited growth of a crop. Huang et al. [11] and Zhao, Chen and Shen [12] empirically assimilated satellite-based remote sensing information into the WOFORST crop model to estimate regional wheat yield. The former used images combined from Moderate Resolution Imaging Spectroradiometer (MODIS) data and three Landsat TM images while the latter used MODIS images only. Launay and Guerif [13] also reported an integration scheme of SPOT satellite images into the SUCROS model to advance its performance for spatial application in prediction of sugar beet production. Although these approaches quantitatively calibrate the crop model to the actual field conditions for each application, it needs achieving the same inputs required for the crop models. The minimum input requirements for simulation include climate data, soil property data, cultivar-specific genetic coefficients, and field management data [8]. These requirements can be too much to allow the model to be run for an appropriate time in some cases.

GRAMI [3,14,15] is a crop model that uses remotely sensed data and is designed to simulate gramineous crops, such as wheat, corn, and sorghum, using simple inputs (i.e., weather and remote sensing data). GRAMI includes a within-season calibration method, which allows the model to fit measured LAI values using an iterative numerical procedure. Thus, model parameters and initial conditions can be adjusted based on comparisons between measured and simulated values. The resulting simulated crop growth minimizes the error between simulated leaf areas and values of leaf areas obtained from remote sensing. An advantage of this procedure is the capability to use infrequent observations to calibrate the model. The GRAMI model was further developed and was applied to simulate cotton growth and lint yield under limited irrigation conditions [16,17], as well as geospatial projections of rice productivity [18]. The crop modelling technique formulated in GRAMI was applied to assess and monitor crop conditions and yields at regional scales, using imagery from operational satellites [19,20,21,22,23]. This within-season calibration methodology was also used to estimate evaporation and biomass production [24,25].

The remote sensing data either from airborne or satellite platforms is spatiotemporal data consisting of many pixels but a relatively small number of time points. For example, the reflectance can be obtained from the RapidEye satellite with a ground resolution of 5 m at 800 × 800 pixels five times during a crop-growing season in a monsoonal climate region. A hybrid approach utilizing both remote sensing and crop modelling can fill the temporal gap. Parameters that characterize a crop model are not necessarily constant. However, pixel-by-pixel estimation is not a good choice for studying the spatial variation in the parameters for the following two reasons. First, the model can be unidentifiable when the sample size is smaller than the number of unknown parameters. Second, as illustrated in James and Stein [26] and Efron and Morris [27,28], the pixel-by-pixel estimate is not optimal with respect to the mean-squared loss. An empirical Bayesian method is introduced to overcome the above-mentioned difficulties. Bayesian methods are characterized by the following concepts and procedures. Random variables are used to model all sources of uncertainty in statistical models, e.g., uncertainty due to lack of information. Determination of the prior probability distribution is required to take into account the prior information. Use of Bayes’ formula is sequentially performed, calculating the posterior distribution when more data become available. In Bayesian statistics, the probability can range from 0 to 1. In this paper, a random-effect model is used to describe the location-dependent unknown parameters, and the empirical Bayesian approach is adopted. A random effect is a factor whose levels are considered a random sample from some population, technically considered to have a normal distribution. Moreover, the relationship between reflectance and LAI is incorporated into the statistical model. These models allow us to monitor the crop growth and estimate the crop model using the remote sensing data from various platforms and optical sensors with different ground resolutions. We report a statistical hybrid approach of crop modelling and remote sensing as well as a method to project spatiotemporal crop growth information using paddy rice (Oryza sativa). It is one of the main food crops for more than half of the world’s population, including about 557 million people in Asia.

2. Materials and Methods

2.1. Study Sites

Ground-measured growth and yield and remote sensing data of crop growth were obtained for three separate rice paddy sites: one (field area of ~0.4 km²) site was located at the experimental fields of Chonnam National University (CNU), Gwangju and the other two sites were commercial farm fields (each field area of ~1 km² at Buan and Gimje), Chonbuk Province, South Korea (Figure 1). The mean annual air temperature and the annual precipitation averaged over the past three decades (1981–2010) are 13.8 °C and 1391 mm year^–1 in Gwangju, 13.3 °C and 1313 mm year^–1 in Gimje, and 12.6 °C and 1250 mm year^–1 in Buan, respectively (Korea Meteorological Administration) [29]. The East Asian monsoon climate is dominant from June to October in these regions, with more than half of the annual precipitation falling during this period.

2.2. Crop Growth Data

Paddy rice growth simulation at the field scale was performed using ground-based remote sensing input data obtained at the CNU site in 2011, 2012, and 2013 and the Buan and Gimje sites in 2014. At the CNU site, 30-day-old rice seedlings (cv. Hopum in 2011, cv. Hwasunchal in 2012, and cv. Unkwang in 2013) were mechanically transplanted in paddy fields on 6 June 2011 (day of year (DOY) 157), 28 May 2012 (DOY 148), and 20 May 2013 (DOY 140). At the Buan and Gimje sites, 30-day-old rice seedlings (cv. Ilmi at Buan and cv. Saenuri and Sindongin at Gimje) were transplanted on 28 May 2014 (DOY 148). Nitrogen (N), phosphate (P), and potassium (K) fertilizers were applied in amounts corresponding to 90 N, 50 P, and 60 K kg·ha^–1, respectively.

Data collected were canopy reflectance of plants, using an MSR16R multispectral radiometer (CROPSCAN Inc., Rochester, Minnesota, USA), and LAI using an LI-2200 Plant Canopy Analyser (LI-COR Inc., Lincoln, Nebraska, USA). These data were measured with four replications on DOY 188, 194, 201, 207, 215, 223, 229, 236, 243, 250, 257, and 265 in 2011; on DOY 173, 184, 193, 201, 207, 213, 226, 237, 243, 249, 262, and 276 in 2012; and on DOY 167, 172, 182, 190, 200, 206, 220, 231, and 241 in 2013 at the CNU site. Data at both the Buan and Gimje sites were measured on DOY 165, 177, 195, 203, 223, 238, 258, and 276 in 2014. The yield was estimated by multiplying four yield components, which were measured three times in the sample plots, based on random sampling at maturity. The yield components were panicle number per m², spikelet number per panicle, percentage of filled grain, and the weight of 1000 grains. Weather data at the CNU site were measured using an automated weather station (WS-GP1, Delta-T Devices, Cambridge, UK). Climate data for simulation for the Buan and Gimje sites were obtained from the Korea meteorological weather stations (https://data.kma.go.kr/cmmn/main.do, accessed on 6 October 2018).

Four vegetation indices (VIs) (

Table 1. Vegetation indices used for rice growth phenological stage monitoring and leaf area index estimation.

Vegetation Index♩	Equation♪
NDVI [30]	$\left(R_{800}-R_{660}\right)/\left(R_{800}+R_{660}\right)$
RDVI [31]	$\left(R_{800}-R_{660}\right)/\sqrt{R_{800}+R_{660}}$
OSAVI [32]	$\left(R_{800}-R_{660}\right)/\left(R_{800}+R_{560}+0.16\right)$
MTVI1 [33]	$1.2\left(1.2\left[R_{800}-R_{560}\left]-2.5\right[R_{660}-R_{560}\right]\right)$

^♩ NDVI = Normalized Difference Vegetation Index, RDVI = Re-normalized Difference Vegetation Index, OSAVI = Optimized Soil-Adjusted Vegetation Index, and MTVI = Modified Triangular Vegetation Index. ^♪ R₈₀₀, R₆₆₀, and R₅₆₀ represent the reflectance values of each waveband (800 nm, 660 nm, and 560 nm).

) were calculated using the reflectance data as input for evaluation of the crop model formulated in this study. The four VIs were empirically chosen from structural indices available and effective for the growth determination of plant canopies using wavebands (i.e., 560, 660, and 800 nm) of multispectral optical sensors in most satellites. In these datasets, both VIs and LAI are observable. However, data of only one point is available.

2.3. Remote Sensing Data

We used remote sensing images from two satellites, i.e., the COMS with the medium ground resolution of 500 m and the RapidEye with the high ground resolution of 6.5 m. The COMS developed by the Korea Aerospace Research Institute (KARI) and launched on 27 June 2010, is a geostationary satellite stationed at an altitude of 3600 km above the Earth’s equator and at a longitude of 128.2°E. The COMS loads two typical sensors of the GOCI and the MI (

Table 2. Detailed characteristics of the Geostationary Ocean Colour Imager (GOCI) and the Meteorological Imager (MI) sensors of the Communication Ocean and Meteorological Satellite (COMS) as well as the Multispectral Push-broom Imager (MPI) of the RapidEye satellite.

Satellite	Sensor and Orbit Type (Altitude)	Channel	Wavelength (µm)	GSD (km)	IFOV (μrad)
COMS	GOCI, Geo-synchronous (36,000 km)	B1	0.40–0.42	0.5	28
		B2	0.43–0.45	0.5	28
		B3	0.48–0.50	0.5	28
		B4	0.55–0.57	0.5	28
		B5	0.65–0.67	0.5	28
		B6	0.68–0.69	0.5	28
		B7	0.74–0.76	0.5	28
		B8	0.85–0.89	0.5	28
	MI, Geo-synchronous (36,000 km)	VIS	0.55–0.80	1	28
		SWIR	3.50–4.00	4	112
		WV	6.50–7.00	4	112
		IR1	10.30–11.30	4	112
		IR2	11.50–12.50	4	112
RapidEye	MPI, sun-synchronous (630 km)	Blue	0.44–0.51	6.5 × 10–3	20
		Green	0.52–0.59	6.5 × 10–3	20
		Red	0.63–0.685	6.5 × 10–3	20
		Red edge	0.69–0.73	6.5 × 10–3	20
		NIR	0.76–0.85	6.5 × 10–3	20

B1 to B6 = visible, B7 and B8 = near infrared, VIS = visible, SWIR = shortwave infrared, WV = water vapor, IR1 = infrared 1, IR2 = infrared 2, and NIR = near-infrared. GSD = ground sampling distance. IFOV = instantaneous field of view.

). The GOCI can observe the Korean landscape eight times a day, with eight spectral wavebands mounted on the system. The primary objectives of the GOCI sensor are: (1) to detect, monitor and forecast short-term biophysical phenomena; (2) to analyse the bio-geochemical variables and the cycles; and (3) to identify information on the yellow dust and land classification. The GOCI was predominantly considered to monitor the ocean phenomena, with a Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) spectral band [34]. However, its high temporal resolution of the observations and vegetation-sensitive spectral wavebands are mostly suitable for land surface-based applications and, more specifically, for monitoring various types of crop information and conditions [23]. In this study, the GOCI is used to determine the cumulative vegetation indices (VIs) in South Korea, since the setup is responsive to vegetation bands (i.e., red and near infrared, NIR).

RapidEye [35] collects 4 million square kilometres of data per day (Table 2). RapidEye satellite images comprise five spectral wavebands (red, green, blue, red edge, and near-infrared) and are provided commercially using three processing levels: “Level 1B” geometrically uncorrected images, “Level 3A” ortho-rectified tile images with radiometric, geometric, and terrain corrections, and “Level 3B” ortho-rectified, bundle-adjusted images that are larger than the Level 3A products [35]. The Gimje and Buan sites were selected, along with the availability of the RapidEye satellite images (Plant Labs, Inc., CA, USA) for 2014, to perform simulations at the regional scale. For the present study, we obtained radiometrically and geometrically corrected Level 3A images, including the Buan and Gimje sites that were acquired on DOY 149 (29 May), 198 (17 July), 221 (9 August), and 256 (13 September) in 2014. These images contained reflectance values with four different wavebands at 11,700 × 7900 pixels with a size of 5 m. The image data were further processed to classify rice paddy field areas using a digitized paddy coverage map from the Korea Ministry of Environment.

Likewise in the model evaluation procedure, four different VIs (Table 1) were determined and used to project spatiotemporal productivity of paddy rice using the developed crop modelling design discussed in the following section.

2.4. Climate Data

The incident solar radiance on the surface (insolation) estimated using the COMS MI and the air temperatures from the Korea local analysis and prediction system (KLAPS) were used to evaluate the updated GRAMI model to simulate crop growth. The insolation data reflected the energy source of photosynthesis in crop canopies. The missions of the MI sensor (Table 2) comprises continuous monitoring of imagery and extracting of meteorological products to allow early detection of severe weather phenomena and to monitor climate change and the atmospheric environment. The MI has been used to estimate the insolation as the infrared channels can be useful for interpreting the complicated cloud effects. In this study, a pixel-based physical model, using instantaneous satellite observations and atmospheric information, was adjusted out of various satellite-based insolation models [36,37,38].

The satellite-based solar insolation from COMS MI was calculated using the Kawamura physical model [39]. This model contains an improved cloud factor, as it considers the visible satellite reflectance and the solar zenith angle instead of the illumination temperature because the passing depth of the cloud is more sensitive to the amount of irradiance attenuation [40]. The physical model used to estimate the insolation is as follows [39,40,41,42]:

$$S_T=S_I+S_R+S_A$$

(1)

where S_T, S_I, S_R, and S_A are the total insolation, direct irradiance, diffuse irradiance due to Rayleigh scattering, and the diffuse irradiance due to scattering by aerosols, respectively. These parameters are used for the physical model of the satellite-based solar insolation [43]. The S_T, S_I, S_R, and S_A formulations are as follows:

$$S_I=S({\tau }_O\cdot {\tau }_R-{\alpha }_w){\tau }_A$$

(2)

$$S_R=S\cdot {\tau }_O(0.5(1-{\tau }_R)){\tau }_A$$

(3)

$$S_A=S({\tau }_O\cdot {\tau }_R-{\alpha }_w)F_C\cdot {\omega }_O(1-{\tau }_A)$$

(4)

$$S=I{(d_M\cdot d^{-1})}^2\cos \theta$$

(5)

where the symbols used are represented as follows: incident solar constant, S; transmittance due to absorption by ozone, τ_o [44]; transmittance due to Rayleigh scattering, τ_R [45]; transmittance due to attenuation by aerosols, τ_A [46]; absorption of water vapor, α_w [44]; ratio of forward to total scattering by aerosols, F_c [47]; single scattering albedo, ω_o; solar constant, I [48]; the annual mean Sun-Earth distance, d_M; Sun-Earth distance, d; and solar zenith angle, θ, respectively.

In this study, the topographical corrections were not considered because the selected crop type was mostly cultivated on the flat land surface, and it was difficult to topographically validate corrected insolation, without using applicable ground measurements on an inclined plane. Most of the pyranometers in South Korea were deployed on a horizontal plane, according to the World Meteorological Organization (WMO) criteria (Guide to Meteorological Instruments and Methods of Observation WMO-No. 8).

The KLAPS used to get the air temperatures was designed to forecast weather conditions of the Korean peninsula with a pixel resolution of 5 km for 12 hours, up to 24 times a day [49]. The KLAPS produces reanalysed data, with a comparatively high-resolution of 1.5 km based on its analytical scheme, using all the possible measured weather data from the region of interest. The KLAPS also adapts the data assimilation part of the local analysis and prediction system (LAPS) developed by the US National Oceanic and Atmospheric Administration/Forecast Systems Laboratory (NOAA/FSL). It is classified into both data collection and analysis modules. The analytical process is composed of the surface analysis procedures and three-dimensional wind, temperature, humidity, cloud, precipitation, and soil analysis procedures. The further detailed procedure of the KLAPS can be found in Albers et al. [50].

2.5. Formulation of the GRAMI Model

The four processes (Figure 2) involved in simulating daily rice growth are: (1) interception and absorption of the incident solar radiation by the leaves; (2) calculation of growing degree-days (GDD); (3) production of the new dry mass by the leaf canopy; and (4) determination of the LAI partitioning of the new dry mass. The details of these procedures and related equations are described in Appendix A. In this study, we applied the same initial conditions and parameter values determined in the earlier studies (

Table A1. Constant and parameter values used for the GRAMI-rice model [18].

Symbol	Description	Unit♪	Value
ε	Radiation use efficiency	g MJ−1	3.49
k	Light extinction coefficient	na	0.6
Ls	Specific leaf area	m2 g−1	0.016
Tb	Base temperature	℃	12.0
L0	Leaf area index at transplant	m2 m−2	0.2
pa	The parameter in the leaf allocation function	na	0.325
pb	The parameter in the leaf allocation function	na	0.00125
pc	The parameter in the leaf senescence function	na	0.00125

^♪ 1 µmol mol⁻¹ is equivalent to 1 ppm.

) [18,21,23].

LAI is a three-dimensional concept, while the reflectance of plants to solar radiation is a two-dimensional concept because it is electronically recorded as a two-dimensional data for the canopies of crops or mostly the top surface of the plants. We presumed that a log-log regression model, with a slope approximately 2/3, could describe the relationship between reflectance and LAI. Based on this theory, the correlations between five VIs (MTVI, NDVI, RDVI, and OSAVI) and LAI were formulated using log-log linear regression models. For each VI, labelled l = 1, 2, 3, 4, and 5, respectively, an empirical model was framed as follows:

$$\log \left({\mathrm{VI}}_{\mathrm{t}}\right)={\alpha }_{VI}+{\beta }_{VI}\log \left(LA{I}_t\right)+{ϵ}_t$$

(6)

where ${\alpha }_{VI},{\beta }_{VI}$, and ${ϵ}_t$ ($\sim N\left(0,{\sigma }_{VI}^2\right)$ represent intercept, slope, and error of the linear regression model, respectively.

The evolution of the LAI for each pixel was explained by the GRAMI-rice model, using four parameters $\mathsf{\theta }=$ (L₀, a, b, and c). These parameters were assumed to be generated from the prior distribution $\mathsf{\psi }\sim \mathrm{N}\left(\mathsf{\mu },\mathrm{D}\right)$, where the transformations:

$$\begin{gathered} \mathsf{\psi }=\left({\mathsf{\psi }}_1,{\mathsf{\psi }}_2,{\mathsf{\psi }}_3,{\mathsf{\psi }}_4\right)=\left(\log \frac{a}{1-a},\log \frac{b}{1-b},\log \frac{c}{1-c},\log \frac{L_0}{1-L_0}\right) \\ \mathsf{\theta }=\mathsf{\theta }\left(\mathsf{\psi }\right)=\left(\frac{e^{{\psi }_1}}{1+e^{{\psi }_1}},\frac{e^{{\psi }_2}}{1+e^{{\psi }_2}},\frac{e^{{\psi }_3}}{1+e^{{\psi }_3}},\frac{e^{{\psi }_4}}{1+e^{{\psi }_4}}\right) \end{gathered}$$

(7)

were used to guarantee that all four parameters (L₀, a, b, and c) range between 0 and 1.

We obtained both the regression coefficients $\left({\alpha }_{\ell },{\beta }_{\ell },{\sigma }_{\ell }^2\right),\ell =1, 2, 3, 4, \mathrm{and} 5$ and the hyper-parameters $\left(\mathsf{\mu },\mathrm{D}\right)$ from the data collected in previous studies [18,21]. These included both the VIs and the measured LAI values. The parameter $\mathsf{\mu }$ was specified using the ‘before-calibration’ values (L₀ = 0.2, a = 3.25 × 10^–1, b = 1.25 × 10^–3, and c = 1.25 × 10^–3). Parameter D is a diagonal matrix with all diagonal elements equivalent to 0.5.

The following numerical procedure was adopted to obtain $\mathsf{\theta }$ for each pixel.

• Step 1: For each pixel, set $\mathsf{\mu }$ served as the initial guess of $\mathsf{\psi }$.
• Step 2: Define ${\mathrm{LAI}}_{\mathrm{t}}=\tilde{\mathrm{G}}\left(t;\psi \right)=G\left(t;\theta \left(\psi \right)\right)$ and consider the objective function:

  $${\displaystyle \sum }_{\ell =1}^5\left\{\frac{1}{{\mathsf{\sigma }}_{\ell }^2}{\displaystyle \sum }_{t=1}^n{\left(\log V{I}_{\ell t}-{\alpha }_{\ell }-{\beta }_{\ell }\log \tilde{G}\left(t;\psi \right)\right)}^2\right\}+{\left(\psi -\mu \right)}^{\text{'}}{D}^{-1}\left(\psi -\mu \right)$$

  (8)
• Step 3: Generate the simulated curve for each pixel from the estimated $\mathsf{\psi }$ in Step 2.
• Step 4: Update $\mathsf{\mu },\mathrm{D}$ as the sample means and sample variances of the estimates in Step 2.

In this procedure, the parameter $\mathsf{\psi }$ was estimated by minimizing the above function, and the optimization was performed using the POWELL optimization routine [51] for one-point simulation cases and the Quasi-Newton minimizer [52] for two-dimentional simulation cases.

2.6. Modelling Spatiotemporal Crop Productivity

The GRAMI-rice model was formulated to integrate remotely sensed data, allowing agricultural system modellers to reproduce and monitor potential crop production information [18]. The GRAMI model can receive remote sensing data as an input to execute the ‘within-season’ calibration procedure [14]. In this process, simulated crop canopy growth (LAI or VIs) is compared with corresponding measured values to allow agreement with the measurement with a minimal error based on parameterization of specified parameters. Four different coefficients (L₀, a, b, and c) are employed in the current model to describe growth processes of rice. Parameter values were obtained through a parameterization process using the Bayesian method with a prior distribution chosen according to the estimates from the previous reports. The relationships between five VIs and LAI were framed using the log-log linear regression models as previously described.

A Crop Information Delivery System (CIDS) was designed earlier as an extended version of the GRAMI-rice model [18] to project pixel-based geospatial crop growth and yield, based on the integration with remote sensing images (Figure 3a). The CIDS employs pixel-based remote sensing data and climate data as the system inputs. The CIDS takes climate data, either from a single weather station or multiple weather stations (pixels) depending on the situation. The GRAMI-rice model is then implemented to simulate crop growth in each pixel using both types of input data.

In the CIDS, the GRAMI model simulates a crop in each pixel using both remote sensing data and climate data as inputs [18]. The new model design was evaluated using paddy field datasets obtained in 2012 at Chonnam National University, Gwangju and in 2014 at Gimje and Gyewha, Chonbuk Province, ROK. The CIDS was then applied using RapidEye satellite images obtained from paddy fields of interest in Chonbuk province, ROK in 2014.

2.7. Statistical Evaluation

Several statistical analytical methods were used to evaluate the reliability of the results. For model evaluation, the mean values of simulated crop yield and LAI were compared with the measured values. A two-sample t-test and two statistical indices of root mean square error, RMSE and model efficiency, ME [53] were applied for the analyses using R software version 3.4.1 [54]. The physical meaning of the ME is that it is a normalized statistic that determines the relative magnitude of the residual variance compared to the observed data variance. ME indicates how well the plot of observed versus simulated data fits the 1:1 line. ME values range from –∞ to 1; the model is more accurate when the value is closer to 1. When values are close to zero, the model predictions are less than or as precise as the observed mean.

3. Results

3.1. Evaluation of Model Performance

The updated GRAMI-rice model performed well in terms of reproducing field conditions of paddy growth in LAI and above-ground dry mass (AGDM). Rice yield could also be determined with reasonable accuracy. For the sake of verification, we used data from CNU, Gwangju from 2011 to 2013 (Figure 4 and

Table 3. Comparisons between simulated and observed mean leaf area index (LAI) of paddy for model calibration at Chonnam National University (CNU), Gwangju and for model validation at Buan and Gimje, Chonbuk in South Korea. Comparison statistics applied were root mean square error (RMSE) and model efficiency (ME).

Division	Year	Cultivar	Simulation	Observation§	RMSE	ME
			------------------ m2 m–2 ------------------
Gwangju	2011	Hopum	3.93	3.84	0.71	0.80
	2012	Hwasunchal	2.98	2.84	0.56	0.86
	2013	Unkwang	2.73	2.69	0.33	0.95
Chonbuk	2014	Sindongjin	4.52	4.57	1.05	0.75
	2014	Saenuri	4.70	4.62	0.99	0.84
	2014	Ilmi	5.04	5.10	0.88	0.86

^§ Observation was based on the GRAMI model projections using four vegetation indices (VIs), i.e., normalized difference vegetation index (NDVI), re-normalized difference vegetation index (RDVI), modified triangular vegetation index (MTVI), and optimized soil adjusted vegetation index (OSAVI).

). Simulated seasonal curves for LAI were fit to the corresponding observed values with a model efficiency (ME) of 0.80 and a root mean square error (RMSE) of 0.71 m² m^–2 for cv. Hopum in 2011, an ME of 0.86 and an RMSE of 0.56 m² m^–2 for cv. Hwasunchal in 2012, and an ME of 0.95 and an RMSE of 0.33 m² m^–2 for cv. Unkwang in 2013. Simulated yields agreed well with the observed yields, with an RMSE of 434.9 kg ha^–1. According to a two-sample t-test (ɑ = 0.05), the simulated yield (μ = 7,365.7 kg ha^–1) was not significantly different (p = 0.883) from the observed yield (μ = 7,190.0 kg ha^–1). Validation, using data from Gimje, Chonbuk in 2014 (Figure 5 and Table 3), showed that simulated LAI values agreed with the observed LAI values with an ME of 0.75 and an RMSE of 1.05 m² m^–2 for cv. Sindongjin, an ME of 0.83 and an RMSE of 0.99 m² m^–2 for cv. Saenuri, and an ME of 0.86 and an RMSE of 0.88 m² m^–2 for cv. Ilmi. Simulated yields also agreed with the observed yields, with an RMSE of 426.5 kg ha^–1. According to a two-sample t-test (ɑ = 0.05), there was no significant difference (p = 0.839) between the simulated yield (μ = 6937.7 kg ha^–1) and the observed yield (μ = 6801.0 kg ha^–1).

3.2. Geospatial Projections of Crop Yield and Growth

The CIDS well projected geospatial variations in rice yield and growth using satellite images from both the COMS with a medium ground resolution of 500 m (Figure 6) and the RapidEye with a high ground resolution of 5 m (Figure 7). When comparing simulated and observed rice yields (Figure 6b), the simulated mean value (μ = 6.77 t ha^–1) was not significantly different (p = 0.392) from the observed mean value (μ = 6.86 t ha^–1), according to a two-sample t-test (ɑ = 0.05). These agreed with an RMSE of 0.44 t ha^–1 and an ME of 0.24. The CIDS was also designed to project crop productivity information with seasonal patterns, as well as two-dimensional variations on any given day of interest during the crop season. The CIDS produced two-dimensional predictive maps of rice grain yield and growth variables well, with a 5 × 5 m pixel resolution (Figure 8). Predicted grain yields ranged between 4150.6 and 8483.2 kg ha⁻¹, with a mean of 6029.1 kg ha⁻¹. The predicted average value agreed with the field measurement of 6381.9 kg ha⁻¹. Spatiotemporal projections of LAI and AGDM reproduced spatial variations of each growth variable well in the scene during the crop season.

4. Discussion

Our study results presented that the GRAMI-rice model was well integrated with either ground-based remote sensing information or satellite-based remote sensing images with different ground resolutions. The modelling regime reproduced crop growth conditions and yields in significant agreement with the corresponding observations. The current study was dedicated to the statistical hybrid approach of crop modelling and remote sensing, as well as the method to project spatiotemporal crop growth information using remote sensing data from various platforms and optical sensors with different geospatial resolutions. On the other hand, the previous studies using GRAMI-rice were focused on applying for different aspects of crop monitoring, most likely as case studies using images from a specific platform, e.g., an unmanned aerial system, UAS [55] or an optical satellite with either a high ground resolution [21] or a coarse ground resolution [23,56,57].

These simulation results demonstrate that the within-season calibration procedure worked well with the observed LAI inputs incorporated into the model. The updated GRAMI-rice model is theoretically incorporated with a Bayesian method for parameter estimation to facilitate an agreement between simulations and observations based on the POWELL [50] or Quasi-Newton [51] optimisation procedures. The POWELL optimisation routine is carried out for one-point simulation cases while the Quasi-Newton minimiser is performed for two-dimensional simulation cases. We assume that the current parameter estimation method is not only technically unique but also advantageous to assimilate various remote sensing data into crop models that have a strong dependence on input LAI from remotely sensed information. The strong integration of the present GRAMI model with remote sensing data can be a discrete advantage in several ways. First, the present approach requires a simple input requirement, which employs only the existing observations that characterize the environmental circumstances. Earlier versions of the current model also showed the corresponding information [14,16,18]. Second, the optimization method allows the model to advance the simulation performance. Third, the GRAMI model can be assimilated with remotely sensed information from various platforms, e.g., a UAS [55] and operational optical satellites on-board different ground resolution sensors [21,56,57]. Finally, once the model is applied to the satellite-based remote sensing, it is applicable for any region of interest on the Earth’s surface. However, limitations include the inadequate representation of remotely-sensed information as well as partial observations available during the crop-growing season. These can ultimately bring about some level of disagreement between simulations and observations.

While some limitations in the GRAMI model exist including a firm reliance on remotely-sensed data needed to achieve the modelling mentioned above, the requirement of input parameters and variables has significant implications, particularly for inaccessible and data-sparse regions. In such regions, this kind of crop models is practically relevant, since it is almost impossible to monitor or reproduce the crop productivity without using operational satellite-based remote sensing information.

The CIDS is designed to project geospatial information of crop growth and yield using the GRAMI model integrated with remote sensing images from various platforms of remote sensing from an UAS [18,55] to various operational optical satellites with different ground resolutions [21,56,57]. There have been similar practical efforts to assimilate a crop model with satellite images to improve the predictive performance of crop yields [11,12,13]. Meanwhile, the CIDS is unique, in that (1) it is formulated as a whole integrated program to project spatiotemporal variations in crop productivity, as well as (2) the CIDS requires simple input parameters and environmental variables utilizing remote sensing images. One of the critical issues to achieve advanced monitoring of crop productivity information using the system would be to determine a consistent representation of the information on crop canopy reflectance. In this regard, one should acquire stable predefined relationships between LAI and canopy reflectance values or VIs of interest for each crop, based on field experiments. It is essential to determine a reliable classification of crops of interest and to establish spatial variations in different planting dates for each crop well to project practical productivity, especially using coarse resolution satellite images. Dealing with these issues is beyond the scope of this study. It is also essential to extract the endmembers of the mixed pixels in the case of using the coarse resolution satellite images. In the case of the fine resolution satellite images, sparse changes of obtained images can be an issue for simulation, using the previous version of the GRAMI model [18]. The earlier GRAMI model required receipt of an even distribution of several input images during the crop-growing season to achieve dependable simulation outputs. This approach has been optimized in the current version of GRAMI based on the advanced empirical approach, using a predefined relationship between LAI and canopy reflectance values of wavebands or VIs of interest.

5. Conclusion

We present functional coupling of crop modelling and remote sensing using an updated GRAMI-rice model that uses remote sensing data, and a CIDS formulated to simulate geospatial rice growth and yield by adapting the GRAMI-rice model for use in this study. The simulated values of rice growth obtained with the parameterized GRAMI-rice model showed good agreement with the corresponding field measurements. We also presented that GRAMI-rice was successfully incorporated with optical satellite data with different geospatial resolutions. Simulated geographical variations in yield were in reasonable agreement with the corresponding observed variations in yield. Therefore, the current study demonstrates that the GRAMI-rice model can be applied to reproducing rice growth and development conditions and productivity based on the integration with remote sensing data from various platforms such as a UAS and different optical satellites with different ground resolutions. The CIDS was applied to monitoring and mapping of rice growth and yield. The GRAMI-rice model has relatively simple environmental input requirements that can be provided by remote sensing. We assume that the CIDS can be potentially applied to crop growth monitoring and yield mapping efforts for croplands of various geospatial scales, ranging from farm fields to regions of interest because remote sensing data can be obtained from observations of small-size unmanned airborne platforms (drones) and existing satellite sensors.