Forecasting Transplanted Rice Yield at the Farm Scale Using Moderate-Resolution Satellite Imagery and the AquaCrop Model: A Case Study of a Rice Seed Production Community in Thailand

Abstract

Thailand has recently introduced agricultural policies to promote large-scale rice farming through supporting and integrating small-scale farmers. However, achieving these policies requires agricultural tools that can assist farmers in rice farming planning and management. Crop models, along with remote sensing technologies, can be useful for farmers and field managers in this regard. In this study, we used the AquaCrop model along with moderate-resolution satellite images (30 m) to simulate the rice yield for small-scale farmers. We conducted field surveys on rice characteristics in order to calibrate the crop model parameters. Data on rice crop, leaf area index (LAI), canopy cover (CC) and agricultural practices were used to calibrate the model. In addition, the optimal rice constant value for conversion of CC was investigated. HJ-1A/B satellite images were used to calculate the CC value, which was then used to simulate yield. The validated results were applied to 126 sample pixels within transplanted rice fields, which were extracted from satellite imagery of activated rice plots using equivalent transplanting methods to the study area. The rice yield simulated using the AquaCrop model and assimilated with the results of HJ-1A/B, produced a satisfactory outcome when implemented into the validated rice plots, with RMSE = 0.18 t ha⁻¹ and R² = 0.88. These results suggest that integration of moderate-resolution satellite imagery and this crop model are useful tools for assisting rice farmers and field managers in their planning and management.

1. Introduction

Rice farming is one of the most critical activities in the global agriculture sector, producing a staple food for the majority of the world’s increasing population, which is expected to exceed 9 billion by 2050 and will require approximately 60% more food. Accurate and up-to-date assessments of the spatial distribution of rice cultivation areas are a vital requirement for all stakeholders including policy makers, rice farmers and consumers. Timely assessment with high precision is also crucial for water resource management, market price control and handling of humanitarian food crises [1,2,3].

To maintain stable rice production, Thailand has recently introduced community-based policies to create regional agricultural plantations. The concept of these policies is to encourage aggregation of small rice fields (farmer-based) into one large (community-based) field. There are three major requirements of this policy: expert field managers, an action plan and site-specific technology [4,5].

In Thailand, where rice farming is commonly practiced in a water-shortage environment, development of management strategies is essential for improving crop production. Under limited water resource conditions, crop models are useful tools for supporting decision-making regarding effective field management. However, such models are complicated and require a large number of parameters. Models also require advanced skills from end-users for model calibration and operation [6,7]. For example, the DASST [8] and APSIM models require many parameters [6,9,10] and the CropSyst, WOFOST and SPAC models are relatively complicated for end-users. These difficulties prevent the development and extensive application of these models [6].

Traditionally, crop monitoring and yield estimation are based on ground surveys and reports but this approach is expensive and time-consuming [11,12,13]. As technologies have developed, crop modeling has become increasingly important for providing the necessary information to farmers. Forecasting crop yield is essential for decision-making regarding soil and plant management, crop selection, marketing, storage and transport [14]. Specifically, forecasting the crop yield is vital for assisting in agricultural planning and decision-making policies. For example, small-scale farmers require information on crop yield to allow them to estimate costs and potential profits before investing their time and money in crop cultivation [15].

Farm-scale yield data have become a primary factor in agriculture policy and crop insurance ratings. Nevertheless, researchers often have trouble analyzing these programs because existing farm-scale yield data are sparse, not necessarily broadly representative and might suffer from selection bias [16]. Crop growth models are useful for estimating growth, development and yield but they need high-quality spatially-distributed input data to derive reliable estimates [17]. Comparison of the simulated rice yield of small-holder farmers and their different agricultural behaviors would provide necessary information for implementing Thailand’s community-based policy. Agricultural monitoring and management have dramatically benefited from the increased availability of a wide variety of remotely-sensed satellite imagery, ground-sensed data (e.g., weather station networks) and crop models, delivering a wealth of actionable data to stakeholders allowing the streamlining and improvement of agricultural practices [18].

The AquaCrop model is another decision-support tool used for modeling and devising strategies for the efficient management of crop-water productivity [19]. This model is also capable of predicting crop productivity, water requirements and water-use efficiency under water-limiting conditions [20]. Moreover, AquaCrop estimated the biological yield of rice more accurately compared to the Crop Environment Resource Synthesis-Rice and ORYZA2000 models [21]. Other previous studies have used the AquaCrop model to simulate crop growth and yield in different environments [22], ranging from potato [23], soybean [24], maize [9,21,22,25,26], wheat [6,20,27,28,29] to rice [30,31,32]. Other studies, however, suggest that AquaCrop be further tested in different locations in order to increase the accuracy of crop production prediction in various conditions [19].

In addition to crop models, remote sensing (RS) is another appropriate tool used to assess crop growth and production at various scales [33,34]. RS data with coarse and middle resolution are widely used in rice cultivation research [35]. Recently, several RS methods have been developed to map rice areas in different parts of the world, such as optical RS-based mapping, microwave RS-based mapping and the integration of the two methods [7]. Using different techniques, optical images from such satellites as the Advanced Very High-Resolution Radiometer (AVHRR) with 1.1 km spatial resolution, the Moderate Resolution Imaging Spectroradiometer (MODIS) with 250, 500, 1000 m spatial resolution and the Landsat-Vegetation data set (30 m spatial resolution) have been used to estimate cropping area of rice and crop yield [20,36,37,38,39,40]. Although these techniques demonstrate advantages for rice monitoring at regional to global scales because of their wide coverage and relatively long data archiving periods, the low resolution of satellite data are not suitable for rice crop mapping in fields that are relatively small in size, irregular in shape and sometimes fragmented by well-built roads and dense water networks, or mostly mixed with other land cover types. These small-scale rice fields are the critical limitation of RS methods because the area of some rice fields is smaller than the spatial resolution of these satellite platforms. As a result, the mixed-pixel problem is prominent, resulting in creating uncertainty in the discrimination of the spectral signatures of rice from other land cover type [40,41]. Few studies have attempted to estimate the rice yield using high-resolution RS data (such as Quickbird; 0.65 m, WorldView; 0.31 m and IKONOS; NIR 3.2 m, PAN 0.82 m) in the past two decades [36] but their approach has encountered problems with swath width and high cost [42,43]. Using Landsat imagery has also encountered a problem with obtaining the cloud-free images due to the temporal resolution, as the revisit frequency is 16 days [7,34]. As a result, it makes it impossible to obtain phenology information during the relevant crop period. Previous studies suggest the combination use of satellite data with an acceptable spatial resolution and a more frequent revisit cycle [40] but again, it would result in high costs.

Despite shortfalls, a recent study found that imagery from the Chinese HJ-1A/B satellite (30 m spatial resolution and 4 days revisiting) has been successfully applied in China to rice-field mapping along with other crops such as rice, maize, sunflower, and wheat [40,44,45,46]. Here in this study, we integrated moderate-resolution optical images from HJ1A/B (30 m) with the AquaCrop model to simulate the yield of transplanted rice on small farms in a Rice Seed Production Community as a pilot case study area to support agricultural planning policies in Thailand. In addition to using Aquacrop model, our study attempted to employ the advantages of HJ1A/B with an adequate spatial resolution for small farm (not less than 900 m² or 1 pixel in the image) and high temporal observation.

2. Materials and Methods

2.1. Study Area and Field Survey

Huaykhamin is a sub-district of Nong Khae district in Saraburi province, central Thailand. The geographic location is 14.42°–14.44° N and 100.87°–100.89° E (Figure 1). Rice is the main crop and most farmers cultivate it twice per year [47].

Thailand is implementing a pilot project aiming to enhance the capacity of rice farmers to reduce planting costs, improve quality and quantity of yields and share knowledge that would allow the implementation of beneficial methods to other paddy areas. This study selected the paddy area of the community of a rice seed production in the Huaykhamin sub-district of Saraburi province, Thailand, as the first pilot of this smart-farm model. Moreover, this area is one of several pilot projects begun in community-based rice farms in recent years. This community has 86 households with plantations of approximately 480 ha in total [48]. Two types of rice cultivation methods are used; namely direct (354 ha) and transplantation (26 ha). This study focused on farmers using transplantation because this method produces the finest seed stock and because the transplanted rice area has occupied more than 40% of paddy area in 2009–2015 of Thailand [49].

2.2. Data Collection and Processing

Two major processes were considered namely the calibration and validation methods. The transplanted rice samples plots were operated as two sets; Group A (23 ha) for the rice plots using the machine and Group B (3 ha) for the rice plots that were manually transplanted. Rice plot vector boundaries were created using the geo-reference Google Earth images to fit with a geometric corrected HJ1A/B image and the topographic map (1:50,000) of the Royal Thai Survey Department(RTSD) [50,51]. At least 20 Ground Control Points (GCPs) were selected on the image and topographic map. Image rectification was based on a first-order polynomial transformation with RMSE of about 0.2 pixel. A flow chart showing the methodology used in this study is presented in Figure 2.

2.2.1. Calibration

The calibrated rice plots were selected for field measurement; those plots were occupied by the same owners on each plantation methods (i.e., machine and manual) as the same agriculture activities, planted date and harvested date. In addition, data from interviews with farmers including rice age, start plantation/expected harvesting date and other agricultural practices were also used for calibration of crop parameter.

Totally, training samples from 8 and 7 plots in the Group A and B respectively were collected (total 35 pixels in the satellite imagery). The field experiment was conducted on 18 sample points for full detailed of field measurement in a selected plot in each group. The field data were collected about weekly in the fields from 13 July to 27 November 2013. The location of rice plot was obtained by a Global Positioning System (GPS) handset. The following data were collected during the field observations; transplanting dates, rice varieties, rice height, leaf area index (LAI referring to the total one-sided leaf area per unit ground area), harvest dates and a number of rice plants. On the harvesting date, 1 m × 1 m sampling plots were measured within the boundaries of both groups to weight the rice and seeds (g). The overall processes of experiment and data collection are shown in Figure 2. Other calibrated plots (same farmer) were included to investigate NDVI and CC relationship (Section 2.5.4) to generate models for NDVI-CC conversion in order to apply with rice plots group of validation (different farmers).

2.2.2. Validation

This process has worked with others 25 plots (126 pixels) for both rice groups, the set of data was collected from different owners. Therefore, the agriculture practice was verities such as date of plantation, harvesting and amount of actual yield. The calibrated crop parameter and formulation of NDVI to canopy coverage were implemented on this sample. The agreeable calibrated crop parameters and NDVI and CC relationship equations were then implemented on these sample pixels as a validated data set.

The total number of pixels which selected for calibration and validation process was excluded ponds, clouds and shadows. There were 161 pixels (approximately 14.49 ha). The calibration of model parameters consisted of 35 training pixels which covering 3.15 ha and 126 test pixels about 11.34 ha for validation.

2.3. Description of the AquaCrop Model

AquaCrop is a multi-crop model developed by the Land and Water Division of the Food and Agriculture Organization (FAO). The model estimates the crop productivity decrease in response to water stress [52]. The AquaCrop model evolved from the proportionality factor between the relative yield decline and the relative reduction in evapotranspiration (Ky) [26] by (1) separating non-productive soil evaporation (E) from productive crop transpiration (Tr); and (2) estimating biomass production directly from actual crop transpiration using a water productivity parameter. These changes lead to the following equation, which is at the core of the AquaCrop model [53]:

B = WP·ΣTr

(1)

where B is the biomass, produced cumulatively (kg per m²), Tr is the crop transpiration (either mm or m³ per unit surface) summed over the period in which the biomass is produced and WP is the water productivity parameter (either kg of biomass per m² and per mm or kg of biomass per m³ of water transpired).

For most crops, only part of the biomass produced is partitioned to the harvested organs to produce yield (Y) and the ratio of yield to biomass is known as the harvest index (HI), hence:

Y = HI·B

(2)

The underlying processes culminating in B and HI are largely distinct from each other. Therefore, separation of Y into B and HI makes it possible to consider the individual effects of environmental conditions and stresses on B and HI separately. The normalized WP (designated as WP*) remains virtually constant over a range of environments. This has fundamental implications for the robustness of the model, which is further enhanced by the daily quantification of HI over the yield formation period. Improved knowledge of plant responses to water stress on short time scales as well as enhanced computational capacity and more accurate procedures to determine daily soil water status, make it possible to simulate yield in daily time steps. This enabled the important transition from a static approach to a dynamic growth model [23,53,54].

It assumes relative evapotranspiration as the key to simulating the yield of a given environment. It divides evapotranspiration into soil evaporation and crop transpiration and integrates a canopy growth and senescence model as the basis for estimating transpiration. The crop yield is computed as a function of the final biomass and the harvest index, while the effects of water stress are isolated to canopy growth, canopy senescence, transpiration and the harvest index [9,54,55].

AquaCrop includes the major components of the soil-plant-atmosphere continuum and the parameters driving phenology, canopy cover, transpiration, biomass production and final yield [54,55]. To run the model, five daily weather input variables are required, including maximum and minimum air temperatures, evapotranspiration, rainfall, annual mean carbon dioxide concentrations of the soil root zone’s water balance and air CO₂ concentration, which affects water productivity and leaf growth [53,55]. The specific features of AquaCrop that distinguish it from other crop models are its emphasis on water and the use of ground canopy cover instead of the LAI, as well as the use of both water productivity values normalized for atmospheric evaporative demand and carbon dioxide concentration, which gives the model an extended extrapolation capacity [54,55]. Compared with other models, AquaCrop is relatively simple to operate by those with little or no research experience. It produces a high level of accuracy and requires only a limited set of input parameters, most of which are rather easy to acquire [20].

2.4. Data Input of AquaCrop Model and Analysis

The first calibration experiment was conducted in both Group A, where the rice cultivar RD49 was mechanically transplanted on 9 July 2013 and harvested on 17 October 2013 and Group B, where the rice cultivar RD31 was manually transplanted on 12 August 2013 and harvested on 22 November 2013. Four types of data were required for the analysis: climate data, soil characteristics, crop characteristics and field and irrigation management. The training points used for calibration of the crop parameters totaled 35 pixels.

2.4.1. Climate Data

A set of climate data such as rainfall, temperature, humidity, wind speed and sunshine were obtained from the local Huaykhamin weather station and input into the daily temporal model. The data were collected during the rainy season (24 June to 31 December 2013). The average rainfall was 1281.6 mm, the average temperature was 27.3 °C, average relative humidity was 79.45%, average wind speed (2 m above the surface) was 0.73 m/s and the average duration of sunshine per day was 6.64 h. The daily reference evapotranspiration (ETo) was computed using the full set of data based on the FAO Penman-Monteith method [56] with the help of ETo calculator software [31,57] and the climate data of the study area as an input. Thus, the maximum reference crop evapotranspiration (ETo) was calculated as 5.1 mm per day, the minimum as 2.3 mm per day and the mean as 3.96 mm per day.

The default CO₂ file supplied with the AquaCrop model was applied for crop yield simulation. Since both the studied groups were located within approximately 500 m of each other, the same climate data were applied to both. Evaporation was simulated using crop transpiration and soil evaporation and the daily transpiration was used to derive the daily biomass productivity of the crop [26].

2.4.2. Soil Data

Entire rice plots of the study area fell within the same soil-series boundary map; Groups A and B were approximately 250 to 300 m from the measurement point located at 14.436° N, 100.878° E as observed by the Land Development Department (LDD) [58]. Soil data were published in 2004 in “Characterization of Established Soil Series in the Central Plain Region of Thailand, Reclassified According to Soil Taxonomy 2003.” The soil series of the study plots is named “Manorom” [59] and the soil horizon profile characteristics are shown in

Table 1. Soil horizon profile characteristics.

Horizon	Description	Thickness (m)	PWP (vol %)	FC (vol %)	SAT (vol %)	TAW (mm/m)	Ksat (mm/Day)	Tau
1	Silty clay	0.16	32.0	50.0	54.0	180	15.0	0.22
2	Silty clay	0.14	32.0	50.0	54.0	180	15.0	0.22
3	Clay	0.30	39.0	54.0	55.0	150	2.0	0.11
4	Clay	0.40	39.0	54.0	55.0	150	2.0	0.11
5	Clay	0.65	39.0	54.0	55.0	150	2.0	0.11

PWP is the permanent wilting point, FC is the field capacity, TAW is the total available soil water in the root zone (mm), Ksat is the saturated hydraulic conductivity (Ksat) and Tau is the dimension drainage characteristic [31]. All values were calculated by the AquaCrop model.

. The Manorom soil series is a member of the fine, mixed, semi-active, isohyperthermic Aeric (Plinthic) Endoaqualfs. They are deep, very strongly acid soils, characterized by a brown loam, silty clay loam, or clay loam A horizon overlying a brown, grading to grayish brown or light brownish gray, argillic B horizon. These soils are mottles throughout with yellowish brown and strong brown coating along root channels in the A horizon and prominent red and yellowish red mottles in the B horizon. These soils have plinthite in the B horizon with few to common spherical hard iron/manganese nodules in the deeper subsoil or throughout the profile [59].

2.4.3. Crop Data

Local crop characteristics such as planting schedule, sowing density, effective rooting depth, flowering and crop germination were collected for model calibration.

To obtain the CC for the input for AquaCrop, we converted the measured LAI to CC because it is the only possible method to obtain CC values during the canopy senescence phase of the crop [53]. During the canopy senescence phase, there is no simple method to measure CC directly, because green and yellow leaves are intermingled and some leaves are partly green and partly yellow.

LAI is one of the most important biophysical parameters characterizing a canopy [60,61,62]. It is a key constraint for carbon assimilation and transpiration rates, which together drive the accumulation of primary crop productivity. Therefore, LAI is commonly used to estimate photosynthesis, evapotranspiration, crop yield and many other physiological processes in agroecosystem studies [63].

The LAI was estimated using the LAI-2000 Plant Canopy Analyzer (LI-COR Inc., Lincoln, NE, USA). LAI measurement began for 9 training sample points (at least 2 times measurement for each point) in both training rice plots of Group A and B for calibration during the tillering stage, when the rice in Groups A and B was 19 and 12 days old, respectively. (Figure 3). The grain yield and rooting depth were determined by removing a 1 m × 1 m frame from each plot. Before the plants were cut, the height of the plants at ground level, the time to emergence, the first day of senescence and the maturity were recorded. The rice plants were randomly measured from 9 points within the plot. The parameters contributing to plant yield were plant height, plant density and the number of rows and lines per square meter.

The AquaCrop model uses a canopy cover (CC) value instead of LAI [6] to calculate the CC during the growing season, the following equation was applied [32]:

CC = 1 − exp^{(−K × LAI)}

(3)

where CC is the canopy cover, K is the extinction coefficient and LAI is the leaf area index. The values of the extinction coefficient for rice were between 0.4 and 0.7 [32]. As the K value has a rather wide range of many rice varieties, one of the purposes of this study was to find an optimal K value for the local rice in the study area. In this experiment, rice extinction coefficients of 0.4, 0.45, 0.5, 0.55, 0.6, 0.65 and 0.7 were investigated before selecting the optimum K value for the transplanted rice parameter calibration.

2.4.4. Irrigation and Field Management

The experimental period was conducted during the rainy season; therefore, we assumed that there was no stress from water usage. Similarly, it was presumed that there was no stress from fertilizer usage because farmers in this community have extensive expertise and knowledge of farm treatments. The soil fertility option was set as non-limiting biomass production, corresponding to soil fertility and the mulch covering the surface was set as 0% according to farmer behavior. However, the field surface practice was set to the soil bund setting, which relates to the volume/area of water storage in the field; this was determined from field observations. Measurement of the bund height was made for all four sides of each rice plot. The average height of the bunds for Groups A and B was 0.49 m and 0.41 m, respectively.

2.4.5. Model Parameters and Calibration

Calibration was performed using conservative crop parameters available in the AquaCrop database. These conservative parameters are crop-specific and do not change materially with time, management practice, geographic location, or climate. Crop parameters, including crop phenology, crop transpiration, biomass production, yield formation and stress, were used to calibrate the model and then to simulate the crop yield. The conservative crop parameter values used in the calibration are listed in the FAO AquaCrop database. Thus, the same conservative parameter values specified for rice in the FAO reference [55,65,66] were input for this study. The other parameters are cultivar-specific or less conservative. The cultivar-specific parameters are affected by location, crop cultivar, the condition of the soil profile and management practices and must be specified by the user.

Table 2. Cultivar-specific parameters input into the AquaCrop model.

Description	Value		Units/Meaning
	Group A	Group B
Soil surface covered by an individual seedling at 90% emergence (cm2/plant)	5	5	CCo
Number of plants per hectare	2,000,000	2,000,000	Plant/ha
Time from transplanting to recover	1	1	days
Maximum canopy cover (CCx)	85.18	90.78	%
Time from transplanting to start of senescence	107	98	days
Time from transplanting to maturity, i.e., length of crop cycle	110	105	days
Time from sowing to flowering	65	65	days
Length of the flowering stage	15	20	days
Minimum effective rooting depth (Zn)	0.30	0.30	m
Maximum effective rooting depth (Zx)	0.35	0.35	m
Reference harvest index (HIo)	36	43	%

describes the cultivar-specific parameters which were observed during the field survey from July to November 2013.

2.4.6. Model Evaluation

Evaluation of model performance is essential for providing a quantitative estimate of its ability to reproduce an observed variable, to evaluate the impact of calibrating model parameters and to compare model results with previous research. Several statistical indicators are available to evaluate the performance of a model. Each has its own strengths and weaknesses, meaning that the use of an ensemble of different indicators is necessary to assess model performance sufficiently.

The AquaCrop model results for both experimental groups were evaluated against the experimental data set from the 2013 rainy season using the calibrated model parameters. To evaluate the goodness of fit between the observed yield and simulated output, statistical indicators such as the coefficient of determination (R²), the Nash-Sutcliffe model efficiency coefficient (EF), Willmott’s index of agreement or compatibility (d), the root mean square error (RMSE) and the normalized root mean square error (NRMSE) were used to compare the simulated and measured values of the parameters [27,67].

R² shows the discrepancy between the simulated and measured values, signifying the proportion of the variance in the measured data described by the model. It ranges from 0 to 1, in which values close to 1 indicate a good agreement and values greater than 0.5 are considered acceptable. EF shows the efficiency of the model in simulations of the parameters. The index of the agreement is a measure of the relative error in the model estimates. It is a dimensionless number and ranges from minus infinity to 1, where 1 indicates a perfect match between the model and the observations. An EF of 0 means that the model predictions are as accurate as the average of the observed data. A negative EF occurs when the mean of the observations is a better prediction than the model [67].

Willmott’s index of agreement (d) is a method used to measure the degree to which the observed data are approached by the predicted data. It ranges between 0 and 1, in which 0 indicates no agreement and 1 indicates perfect agreement. The root mean square error (RMSE) is the most widely used statistical indicator for measuring the average magnitude of the difference between predictions and observations. It ranges from 0 to positive infinity. The normalized root mean square error (NRMSE) is the normalized RMSE using the mean of the observed variable $\overline{\mathrm{O}}$. The NRMSE is expressed as a percentage and gives an indication of the relative difference between the model and the observations. A simulation can be considered excellent if the NRMSE is smaller than 10%, good if between 10% and 20%, fair if between 20% and 30% and poor if larger than 30% [67].

2.5. Remote Sensing Data and Analysis

RS techniques provide an appropriate tool for delineating spatiotemporal patterns of crop status on a per-pixel basis. The spectral characteristics of the data may be affected by a mixed-pixel problem due to the relatively coarse resolution of spatial and temporal data [41]. Therefore, we used moderate-resolution HJ-1A/B CCD imagery for analysis because this satellite constellation can revisit the same location within two days, resulting in higher temporal precision for long-term dynamic rice monitoring [68].

2.5.1. HJ-1A/B CCD Images

The three small sun-synchronous HJ-1 satellites used for environmental and disaster monitoring and forecasting in China were launched in 2008. Both HJ-1A and HJ-1B employ CCD cameras, while the former includes an infrared camera and the latter a hyperspectral camera (HJ-1C employs an S-band SAR package which was not used in this study). These satellites have a spatial resolution of 30 m and a revisit cycle of four days (the revisit cycle of the constellation is 2 days) [40,41]. Among the payloads aboard two satellites, the multispectral CCD cameras are important instruments that are widely used in eco-environment monitoring. Each satellite carries two CCD cameras, named CCD1 and CCD2. The HJ-1A/B CCDs have three visible bands (0.43–0.52 μm, 0.52–0.60 μm and 0.63–0.69 μm) and one near-IR and 0.76–0.90 μm) [61]. This study applied the multispectral CCD camera of visible band 3 and near-IR band 4 for vegetation index analysis (

Table 3. Technical specifications for the HJ-1A/B satellites.

Satellite	Payload	Band No.	Spectral Range (μm)	Spatial Resolution (m)	Swath Width (km)	Repetition Cycle (Day)
HJ-1A/B	Multispectral CCD camera	1	0.43~0.52	30	360 (700 for two)	4
		2	0.52~0.60	30
		3	0.63~0.69	30
		4	0.76~0.90	30

). For this study, HJ-1A/B data were acquired from July to November 2013 for a total of 28 images (

Table 4. Dates of the selected HJ-1A/B imagery and field observations.

No.	Satellite	Sensor	Date	Quality	Field Observation Date
1	HJ-1B	CCD1	2 July 2013	7	-
2	HJ-1B	CCD2	5 July 2013	7	-
3	HJ-1B	CCD1	6 July 2013	4	-
4	HJ-1B	CCD2	13 July 2013	7	23 July 2013
5	HJ-1A	CCD2	26 July 2013	2	27 July 2013
6	HJ-1A	CCD2	3 August 2013	3	5 August 2013
7	HJ-1A	CCD1	17 August 2013	4	19 August 2013
8	HJ-1B	CCD2	28 August 2013	8	26 August 2013
9	HJ-1A	CCD1	3 September 2013	7	2 September 2013
10	HJ-1A	CCD1	11 September 2013	6	9 September 2013
11	HJ-1A	CCD2	26 September 2013	3	23 September 2013
12	HJ-1B	CCD1	2 October 2013	4	30 September 2013
13	HJ-1B	CCD1	6 October 2013	5	7 October 2013
14	HJ-1B	CCD2	9 October 2013	8	9 October 2013
15	HJ-1A	CCD2	11 October 2013	8	14 October 2013
16	HJ-1A	CCD1	27 October 2013	8	27 October 2013
17	HJ-1A	CCD1	31 October 2013	8	-
18	HJ-1B	CCD2	1 November 2013	8	-
19	HJ-1B	CCD1	2 November 2013	8	-
20	HJ-1A	CCD2	3 November 2013	9	-
21	HJ-1B	CCD2	5 November 2013	8	4 November 2013
22	HJ-1B	CCD1	6 November 2013	7	-
23	HJ-1B	CCD1	9 November 2013	7	-
24	HJ-1B	CCD2	13 November 2013	7	-
25	HJ-1A	CCD2	15 November 2013	7	12 November 2013
26	HJ-1A	CCD2	25 November 2013	7	19 November 2013
27	HJ-1A	CCD1	27 November 2013	7	26 November 2013
28	HJ-1A	CCD2	30 November 2013	8	-

) with mostly ≥ quality of 7 (relatively to cloud cover ≤30%). Table 4 also shows the source satellites (1A or 1B), sensors (CCD1 and CCD2), the date of the satellite acquisitions, quality and files observation date. For the Solar elevation and azimuth can be collected from the image header files. The field observation dates, which occurred 17 times during one crop cycle. Even though some full scene of the images obtain lots of clouds more than 50% but as the study area rather small therefore some images were not hidden with cloud and shadow.

2.5.2. Satellite Image Pre-processing

The HJ-1A/B images acquired from July to November 2013 underwent pre-processing. The atmospheric correction consisted of three steps. First, the digital number (DN) values were converted to at-sensor spectral radiance (L_sat) by radiometric calibration. G (Wm⁻²·sr⁻¹·μm⁻¹) is the calibration factor and L_o (Wm⁻²·sr⁻¹·μm⁻¹) is the calibration offset. G and L_o were restored in the corresponding header file and are displayed in

Table 5. Coefficients used for radiometric calibration and mean extraterrestrial solar irradiance (E_sun, Wm²·μm⁻¹) used for TOA reflectance of HJ-1A/B CCD 2013 images.

Satellite Sensor		Band 3		Band 4
		G/Lo	Esun	G/Lo	Esun
HJ-1A-CCD1	G	1.0673	1542.664	1.0429	1073.826
	Lo	3.6123		1.9028
HJ-1A-CCD2	G	1.0464	1549.824	1.0519	1078.317
	Lo	3.7360		0.7385
HJ-1B-CCD1	G	1.0555	1566.714	1.1042	1077.085
	Lo	3.2144		2.5609
HJ-1B-CCD2	G	1.0245	1553.201	1.0146	1074.544
	Lo	5.2537		6.3497

. Second, the resulting images (L_sat) were converted into exo-atmospheric top-of-atmosphere (TOA) reflectance.

Third, the exo-atmospheric TOA reflectance (P_TOA) was converted into surface reflectance (P_λ) using the atmospheric code of the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) model, a radiative transfer model used to process visible and near-infrared multi-angle data [69,70].

After the surface reflectance conversion, the images were geometrically corrected using a WGS-84 projection using the geometric corrected HJ-1A/B based image as it is recommended by the CSRS center [50]. The original of the based image was referenced to the topographic map 1:50,000 of RTSD [51]. The GCPs (at least 20 points over the sub-subset image scene with coverage of about 800 km²) were selected throughout using apparent features such as buildings and crossroads. To minimize the statistical properties of image alteration, nearest neighbor resampling was applied for all geometrical transformations. The polynomial transformation method was used to re-project the images with an RMSE of less than 0.3 pixel.

2.5.3. Vegetation Index Calculation

The surface reflectance normalized difference vegetation index (NDVI) is most commonly applied to optical images for crop information [7]; the resulting values can then be calculated for leaf area index and canopy cover conversion. NDVI is calculated as:

$$\mathrm{NDVI}=\frac{{\mathrm{P}}_{\mathrm{NIR}}-{\mathrm{P}}_{\mathrm{RED}}}{{\mathrm{P}}_{\mathrm{NIR}}+{\mathrm{P}}_{\mathrm{RED}}}$$

(4)

where P_NIR is the surface reflectance in the near infrared band and P_RED represents the surface reflectance in the red spectral band.

2.5.4. Regression Analysis for NDVI and LAI/CC

Empirical regression analysis was performed on the LAI field measurements for both rice groups and the average LAI was used to explore the general pattern of seasonal LAI in the training rice plots. The NDVI to CC conversion was accomplished using an empirical relationship between the NDVI from the HJ-1A/B imagery and the CC (in situ LAI converted to CC using Equation (3)). The regression model was evaluated using trend analysis and the coefficient of determination (R²).

2.6. Validation of Simulated Yield

After interpolation of the CC from the NDVI of all HJ-1A/B images, model results using the best-fit polynomial equations for training points of Groups A and B were applied to 126 tested pixels in 25 rice plots owned by various farmers for validation. However, the pixels were grouped as per the two transplanted methods. The polynomial function was selected because the NDVI can provide information on the dynamics of vegetation growth. From rice sowing to harvesting, the NDVI of crops first increases with vegetation growth, reaches a peak value at peak growth, then gradually decreases with senescence [71].

As a result, the polynomial function can also be used for CC interpolation. The simulated maximum CC (CC_x) was used as a user-specific input of the 126 tested pixels for validation. The calibrated of the cultivar-specific parameters in Table 2 were used as the input crop parameters along with the other conservative parameters (e.g., climate, soil and the model constant values) the for the validation process as well. The input and output dates of the simulation were set according to the actual transplanting and harvest dates, which were established by farmer interviews otherwise by satellite imagery analysis. Finally, the error of yield simulation was evaluated using statistical indicators.

3. Results

3.1. Calibration

3.1.1. Leaf Area Index (LAI) Measurement

The LAI was measured in the study area during rice cultivation in the rainy season of 2013. LAI values from the calibrated plots were measured using the LAI-2000 Plant Canopy Analyzer. (LI-COR Inc., Lincoln, NE, USA). Figure 4 shows a good relationship between LAI and rice age (days) of a primary rice cycle, which fits with the linear function. At the average of LAI value, R² is approximately 0.917 for Group A (Machine transplant) and 0.88 for Group B (Manually transplanted rice).

In additional, using the One-Sample Kolmogorov-Smirnov test of the distribution of LAI data from field measurement was run. As a result, it presented the normal distribution for LAI data of both rice group A (test statistic = 0.167, N = 17, p > 0.05) and B (test statistic = 0.189, N = 17, p > 0.05).

A Pearson’s correlation was run to determine the relation between average LAI values and rice age (growth stages in 17 field-observed dates). There were very strong, positive correlation between LAI and rice age (days) test both rice group A (r = 0.958, N = 17, p < 0.001) and B (r = 0.941, N = 17, p < 0.001).

3.1.2. Canopy Cover (CC)

Equation (3) was applied to convert the LAI value measured from the field to CC. The values of the extinction coefficient for rice are between 0.40 and 0.70 [32]. Therefore, an experiment to select the optimal extinction coefficient (K) value was performed. The different K values were tested for 7 values from 0.4 to 0.7. The average percentage of CC values from field observations was compared with the average percentage of CC values from the AquaCrop model simulation and the results are shown in Figure 5 and Figure 6 and

Table 6. Statistical indicators used for evaluating the CC simulation performance of the model using different K values for Groups A and B.

Rice Farm	K Value	Canopy Cover (CC)			R2	RMSE	NRMSE	EF	d
		Observation/Simulated Sets (n)	Average Observed CC (%)	Average Simulated CC (%)
Group A	0.40	17	42.90	54.40	0.86	13.9	32.40	0.50	0.86
	0.45	17	46.40	57.40	0.88	13.4	28.90	0.59	0.89
	0.50	17	49.20	60.40	0.89	13.5	27.40	0.62	0.90
	0.55	17	52.00	62.70	0.90	13.0	24.90	0.67	0.91
	0.60	17	54.60	64.90	0.91	12.6	23.10	0.71	0.92
	0.65	17	56.80	66.40	0.92	11.9	20.90	0.76	0.94
	0.70	17	58.90	67.80	0.93	11.2	19.10	0.79	0.94
Group B	0.40	17	51.10	54.90	0.88	9.10	17.80	0.85	0.96
	0.45	17	54.40	57.50	0.89	8.7	16.00	0.88	0.97
	0.50	17	57.30	59.50	0.90	8.3	14.50	0.89	0.97
	0.55	17	59.30	61.50	0.92	7.7	13.00	0.91	0.98
	0.60	17	62.20	62.80	0.92	7.8	12.50	0.92	0.98
	0.65	17	64.30	64.20	0.92	7.6	11.90	0.92	0.98
	0.70	17	66.10	65.50	0.93	7.6	11.50	0.92	0.98

. For both groups, the best model performance was derived with a K value of 0.70.

3.1.3. Evaluation of the Calibrated Crop Parameters

The AquaCrop model was calibrated using parameters that concentrated on the local crop characteristics, including the conservative and user-specific parameters. The details of user-specific parameters are shown in Table 2. The model performance was evaluated using the statistical indicators shown in

Table 7. Statistical Results of Yields Simulated from the Calibrated Crop Parameters of the Training Points.

Sample Rice Farm	Yield (t ha−1)			EF	d	RMSE	NRMSE
	Observed (Seed Sample Weighting)	Observed (Average/Interview)	Simulated (AquaCrop)
A	7.72	6.00	6.39	0.40	0.60	1.33	22.15
B	7.17	5.00	7.53	0.97	0.99	0.36	7.22

, which explain the yield values observed from the calibrated training points of the two groups. The observed yield of a 1 × 1 m (seed weighting) in each rice sample plots and the average observed yield (determined by interview) for Group A was 7.72 and 6.00 ton/ha, respectively and 7.17 ton/ha and 5.00 ton/ha for Group B, respectively. In addition, each farmer uses private service of a rice combine harvester for all of their own rice plots. Then the rice weighting process at a rice mill. Therefore, farmers get the rice weight recording then they average as per their occupied rice area; kg per Rai (kg per 0.16 ha). The yields simulated by the AquaCrop model for Groups A and B were approximately 6.39 ton/ha and 7.53 ton/ha, respectively.

3.1.4. Relationship between Surface Reflectance NDVI and CC

The results of the regression analysis between surface reflectance NDVI (from the HJ-1A/B images) and CC establishment (using Equation (3)) for Groups A and B are shown in Figure 7. The selected model equations for the regression analysis are summarized in

Table 8. Summary of testing points (126 pixels) used for validation of the AquaCrop model’s calibrated crop parameters.

Rice Group	Plot ID (Plots)	Plot Area (ha)	Number of Validation Points (Pixels)	Test Pixels Area (ha)	Transplant-Harvest Date	Avg. CCx (%)	Observed Yield (t ha−1) (Interview)	Simulated Yield (t ha−1)	ΔYield (t ha−1)
A	1a	0.77	4	0.36	9 July–6 November	79.05	6.00	6.21	0.21
	1b	0.38	3	0.27	9 July–27 October	79.09	6.00	6.13	0.13
	1c	2.05	11	0.99	9 July–27 October	79.08	6.00	6.13	0.13
	1d	2.22	9	0.81	9 July–27 October	79.06	6.00	6.13	0.13
	2a	1.51	9	0.81	26 July–7 November	79.07	5.00	5.16	0.16
	2b	1.10	6	0.54	26 July–7 November	79.04	5.00	5.16	0.16
	2c	0.38	1	0.09	26 July–7 November	78.68	5.00	5.16	0.16
	3a	0.90	5	0.45	18 August–30 November	79.09	6.25	6.12	0.13
	3b	0.78	3	0.27	15 July–30 November	79.10	6.25	5.83	0.42
	3c	0.89	4	0.36	15 July–30 November	79.08	6.25	5.83	0.42
	3d	0.74	2	0.18	18 August–30 November	79.10	6.25	6.12	0.13
	3e	0.67	3	0.27	18 August–30 November	79.10	6.25	6.12	0.13
	3f	0.31	1	0.09	18 August–30 November	79.10	6.25	6.12	0.13
	3g	0.33	1	0.09	15 July–22 November	79.08	6.25	5.83	0.42
	3h	0.36	1	0.09	15 July–22 November	79.09	6.25	5.83	0.42
	4a	1.91	9	0.81	10 July–15 December	79.04	6.19	6.13	0.06
	5a	1.09	5	0.45	10 August–28 November	79.07	6.00	6.24	0.24
	5b	0.36	2	0.18	10 August–28 November	79.05	6.00	6.24	0.24
	6a	2.12	7	0.63	15 July–6 November	79.05	6.00	6.13	0.13
	6b	1.26	4	0.36	15 July–6 November	79.02	6.00	6.13	0.13
	7b	2.77	17	1.53	20 July–16 November	79.06	6.00	6.19	0.19
Sum. Rice Group A =	21	22.89	107	9.63	Min. =	78.68	5.00	5.16	0.06
					Max. =	79.10	6.25	6.24	0.42
					Avg. =	79.06	5.91	5.98	0.17
					Stdv. =	0.06	0.40	0.36	0.09
					MAE =	-	-	0.175	-
					RMSE =	-	-	0.194	-
B	8a	1.09	8	0.72	10 August–20 November	91.98	5.00	5.10	0.10
	8b	1.17	7	0.63	10 August–12 November	92.24	5.00	5.10	0.09
	8c	0.48	3	0.27	10 August–12 November	92.33	5.00	5.07	0.07
	8e	0.33	1	0.09	10 August–12 November	92.53	5.00	5.10	0.10
Sum. Group B =	4	3.07	19	1.71	Min. =	90.49	5.00	5.01	0.01
					Max. =	92.55	5.00	5.10	0.10
					Avg. =	92.16	5.00	5.07	0.07
					Stdv. =	0.50	0.00	0.02	0.02
					MAE =	-	-	0.071	-
					RMSE =	-	-	0.074	-
Total validation rice all plots	25	25.96	126	11.34	Average =	80.87	5.80	5.87	0.16
					MAE	-	-	0.159	-
					RMSE =	-	-	0.182	-
					R2 =	-	-	0.882	-

. The graphs show the best fit with the polynomial function for both groups, as well as for the relationship between NDVI and %CC. The highest R² values for the polynomial formula of NDVI to %CC for group A (average surface reflectance NDVI from 26 training pixels and B (average surface reflectance NDVI from 9 training pixels) were 0.94 and 0.86, respectively.

Besides, using the linear relationship (Pearson’s correlation) was operated to determine the relationship between the average NDVI and CC which collected for 17 field-observed dates. Nevertheless, as some images have cloud covered on the experiment plots, therefore those pixels were eliminated (11 images were applied here). There were very strong, positive correlation between NDVI and CC of rice group A (r = 0.860, N = 11, p < 0.001) and strong positive correlation between two variables on rice group B (r = 0.713, N = 11, p < 0.001).

3.2. Validation

The polynomial equations for Groups A (Figure 7a) and B (Figure 7b) were selected to generate the CC of 126 testing pixels for validation of the simulated yield. Their maximum canopy cover (CC_x) was calculated and input into the AquaCrop model. The other crop parameters (crop phenology, biomass production and yield information) were implemented using the calibrated crop parameters presented in Table 2. The dates of transplanting and harvesting were determined by interviews. The simulated yields were compared with the average yield (t ha⁻¹) recorded from farmer interviews (farmers also recorded in farm-report of the Provincial Agricultural Office) for entire 25 transplanted rice plots in the study area. The error of the simulated yield (ΔYield) compared to the observed average yield for Group A was approximately 0.17 (t ha⁻¹) and the RMSE was 0.194 t ha⁻¹; for Group B, the error of the simulated yield (ΔYield) was 0.07 t ha⁻¹ and the RMSE was 0.074 t ha⁻¹. The total values for all validated rice plots (Groups A and B) were an average ΔYield of 5.87 t ha⁻¹, RMSE of 0.182 t ha⁻¹ and R² of 0.88. Mean Absolute Error (MAE) of Group A was 0.17 and 0.07 for Group B, while total both groups was 0.16 (Table 8). Moreover, the result of Pearson’s correlation calculated between simulated yield and the observed yield showed a positive correlation of the validated test pixels (r = 0.939, N = 126, p < 0.001).

4. Discussion

RS data contain the growing condition of the crop at the observation time only, making it difficult to estimate crop yield directly [72]. In addition, although RS-based empirical forecasting models are relatively simple to build, these models cannot take into account temporal changes in crop yields without long-term field experiments [7]. Recently, the integration of RS data into crop growth models has been investigated as a way to estimate crop yield, becoming increasingly recognized as a potential tool for yield forecasting [7,72]. Moreover, various types of satellite imagery need to be considered to fit the intended utilization and resolution (e.g., spatial, spectral and temporal). However, the limited synergy between RS-based methods and crop models is complicated and so many input parameters are required, including several biophysical parameters (e.g., soil and meteorological variables) and plant parameters (e.g., biomass, LAI and height, age, etc.) which are usually expensive, labor-intensive and time-consuming to acquire [6,7]. In this study, we combined the advantages of the AquaCrop model [20,40] and moderate-resolution HJ1A/B data to eliminate such limitations, achieving successful results with satisfactory validation.

4.1. Canopy Cover

The canopy cover (CC), detected as LAI from the RS data, is one of the most important variables for crop models [60,61] but is hard to determine in the field. According to field experiment results (Figure 4), there is a good relationship between the measured average LAI and rice age for the two experimental rice groups. This explains why the LAI value increases with rice age/growth stage. The leaf area of rice plants increases with growth, reaches a maximum around heading and decreases thereafter due to the death of lower leaves [73]. The presence of this relationship shows that the field data collection was highly accurate for crop parameter calibration.

When LAI was converted to CC using Equation (3), the calibrated best K value for both experimental rice plots was 0.70, signifying that this is the optimal constant value to use when applying the LAI to the CC conversion formula for similar local rice cultivars. The maximum CC was set as a crop parameter along with other cultivar-specific parameters collected from the field (Table 2). The model was calibrated with this set of parameters, producing satisfactory simulated yields (Section 3.1.3). Consequently, the HJ1A/B images were selected to transform values to CC, using the regression equations (Figure 7) and polynomial formulae with higher R² and graph relationships closer to crop growth phenology [74]; therefore, these were proposed for the surface reflectance NDVI to CC conversion model. With the previously calibrated parameters (35 pixels in 15 plots), the CC values from HJ1A/B data were substituted and the yield was simulated in 126 test pixels in 25 plots for validation (Section 3.2). Furthermore, the maximum CC (CC_x) percentage was more constant in the plots transplanted by machine. The validated rice plots (group A) showed only slightly different CC_x values (average CC_x 79.06%, standard deviation 0.06), while the manually transplanted rice plots (group B) had fluctuated CC_x; the average CC_x was 92.16% and the average standard deviation was 0.5. The LAI and CC were calculated using the surface reflectance NDVI of the HJ-1A/B images, producing significant crop phenology that could be applied in future crop models as well as the image classification for rice cultivation methods.

4.2. Yield Simulation

After the application of calibrated crop parameters for yield simulation, there was a strong relationship between the observed and simulated yields in Group B, which used manual transplanting. The Nash-Sutcliffe model efficiency coefficient (EF) was 0.97, indicating a near-perfect match between the observed and simulated yields. Willmott’s index of agreement or compatibility (d) is close to 1, which signifies an almost-perfect agreement between the observed and simulated yields. The RMSE of rice group B was 0.36, indicating a few errors. The NRMSE is 7.24%, indicating an excellent relationship between the measured and predicted yields. The simulated yield of Group A was evaluated as an acceptable result with statistical indicators of 0.40%, 0.60%, 1.35% and 18.99% for EF, d, RMSE and NRMSE, respectively. Nevertheless, the validation of the simulated yields presented good results (MAE = 0.159, RMSE = 0.182 t ha⁻¹ and R² = 0.88) as shown in Table 8 and Figure 8.

5. Conclusions

In this study, we presented a new approach for rice yield estimation at the farm level by integration of moderate-resolution optical satellite imagery and a crop model. The AquaCrop model was proposed as one of the best models for rice yield prediction as it requires a lesser number of input parameters in comparison to other crop models. The RS data provided useful crop information for parameter calibration and canopy cover (CC) calculation. One benefit of this approach is the reduction in cost of field experiments/observations because the moderate-resolution imagery provides enough quality for crop information on small rice farms. The results show that the performance of the calibrated AquaCrop model was satisfactory compared to observed crop parameters from the small paddy areas, adequately simulating rice yield in the rainy season for individual farms with different planting dates, planting methods and rice varieties. Although small numbers of cultivar-specific crop parameters were applied for calibration, the simulated yield results were good.

The outcome set of the cultivar-specific parameters which were calibrated in this research could be applicable for initial model calibration in other paddy areas. Although the range of the constant value can be scaled down for model calibration before application, it would be necessary to adjust the user-specific parameters based on different regions and conditions. Furthermore, we found that CC calculation using surface reflectance NDVI from HJ-1A/B images, HJ1A/B can deliver useful information for the yield forecasting model and achieved satisfactory results in this case study of rice grown using the transplant method.

We conclude that assimilation of moderate-resolution HJ1A/B imagery can be efficiently used with the AquaCrop model to extract crop information that would otherwise be difficult and costly to collect by field observations. This method can, therefore, be suggested to farmers or field managers for use as part of community-based policy and other agriculture planning projects. We further suggest that this model can be implemented in other paddy areas. To increase even more accuracy of our study, future research should focus on the application of rice crop cycles to the dry season in continuous crop-cycle years, crop parameters (LAI, CC, agricultural practices, etc.) for more cultivar rice varieties and the number of test plots/area in order to improve the simulated yield accuracy and expand applicable areas. Concerning a different kind of the method of the rice plantation would be challenged in the future study. With regards to the pre-processing of HJ1A/B satellite data, we also suggest that the importance of spatial autocorrelation should be taken into account for improved classification accuracy in future studies. The future study could evaluate the optimize of various vegetation index for LAI/CC extraction that would be developed and integrated with radar imagery. Also, cloud eliminate techniques should improve the results on multi-aspect.

In addition, to increase the potential of yield forecasting and implementation, long-term meteorological data should be sufficiently developed in rice growing areas and countries for reliable estimates of rice production.