Remote sensing-based estimation of rice yields using various models: A critical review

ABSTRACT Reliable estimation of region-wide rice yield is vital for food security and agricultural management. Field-scale models have increased our understanding of rice yield and its estimation under theoretical environmental conditions. However, they offer little information on spatial variability effects on farm-scale yield. Remote Sensing (RS) is a useful tool to upscale yield estimates from farm scales to regional levels. Much research used RS with rice models for reliable yield estimation. As several countries start to operationalize rice monitoring systems, it is needed to synthesize current literature to identify knowledge gaps, to improve estimation accuracies, and to optimize processing. This paper critically reviewed significant developments in using geospatial methods, imagery, and quantitative models to estimate rice yield. First, essential characteristics of rice were discussed as detected by optical and radar sensors, band selection, sensor configuration, spatial resolution, mapping methods, and biophysical variables of rice derivable from RS data. Second, various empirical, process-based, and semi-empirical models that used RS data for spatial estimation of yield were critically assessed – discussing how major types of models, RS platforms, data assimilation algorithms, canopy state variables, and RS variables can be integrated for yield estimation. Lastly, to overcome current constraints and to improve accuracies, several possibilities were suggested – adding new modeling modules, using alternative canopy variables, and adopting novel modeling approaches. As rice yields are expected to decrease due to global warming, geospatial rice yield estimation techniques are indispensable tools for climate change assessments. Future studies should focus on resolving the current limitations of estimation by precise delineation of rice cultivars, by incorporating dynamic harvesting indices based on climatic drivers, using innovative modeling approaches with machine learning.


Introduction
The world's population is projected to increase to 9.8 billion by 2050 (UN-DESA 2017), demanding more food production. Rice (Oryza sativa L.) is a significant staple crop for billions of people (Purugganan 2014;Dugan 2015). However, climate change will adversely affect future production of rice Hoekstra and Mekonnen 2012). To meet the needs of the projected population, there is a need to reliably estimate rice cultivation area and yield.
Rice yield is defined as grain weight at 14% moisture per unit area for paddy (rough) rice or with the husk attached (Yoshida 1981). Traditionally, rice yield is estimated by interviewing farmers or via field surveys. Both methods are time-consuming, costly, and subject to bias, with under-and overestimations frequently encountered (Fermont and Benson 2011;Sud et al. 2016). On the other hand, advances in geospatial technology like Remote Sensing (RS) and Geographical Information Systems (GIS) now provide potential alternatives for reliable estimation of rice yield in a cost-effective and timely way over regional scales.
Remote sensing provides useful information on rice, including its seasonal and areal extent. In addition to mapping, remote sensing has been used increasingly to estimate rice cultivation area (Asilo et al. 2014;Kwak et al. 2015;Zhang et al. 2017). New sensors, algorithms, computing power, and crop simulation can all potentially improve the accuracies of yield estimation. Previous literature showed wideranging applications of geospatial data and techniques to estimate yield for a variety of crops (Kuenzer and Knauer 2013;Mosleh, Hassan, and Chowdhury 2015;Dong and Xiao 2016;Jin et al. 2018;Dorigo et al. 2007; Kasampalis et al. 2018). For comprehensive reviews of rice area mapping methods, the reader is directed to Kuenzer and Knauer (2013), Mosleh, Hassan, and Chowdhury (2015), and Dong and Xiao (2016). Despite these exhaustive studies, there is no comprehensive study focusing on rice and geospatial yield estimation techniques. This paper aims to synthesize significant concepts in rice yield estimation using RS and geospatial methods, identify knowledge gaps, and point out the direction of future research in this area. Figure 1 illustrates concepts tackled in this paper. We first provide an overview of current literature about RS characteristics of rice and methods useful for representing rice yield spatially. After critically assessing current methods, this paper makes recommendations about how to overcome these limitations. Ultimately, we ask, "how can spatial rice yield estimates be made more accurate?" This review, thus, aims to contribute the literature about how to operationalize reliable and timely estimation of rice yields by consolidating knowledge and proposing the ways forward.

Remote sensing of paddy fields
Remote sensing of paddies involves detection of biophysical variables and mapping of their areal extent. The spectral response of various surfaces found in rice fields is exploited by sensors to derive vital physiological indicators of vegetation health or development stages-the key to determine final yield, which we will discuss in more detail in Section 2.1. Sensor type and their spatial and temporal resolutions also affect yield estimates, as paddies are diverse ecosystems in terms of size, management and environmental conditions; thus, Section 2.2 details the range of sensors previously used. Yield models focus on the behavior over time of a single or several variables to link remotely sensed information to biophysical parameters, outlined in Section 2.3. For instance, Leaf Area Index (LAI), the quantity of leaf area per unit ground area, can be derived from remotely sensed data. It is a good predictor of yield and incorporated into numerous yield estimation models as a state variable Keating et al. 2003). Lastly, accurate rice area mapping is crucial as yield is an area-dependent quantity, so Section 2.4 discusses new trends in mapping rice area with the advent of new technologies such as GEOBIA, sensor fusion and machine learning, in addition to traditional pixel-based methods.

Life cycle of rice
Understanding rice growth stages is critical to develop reliable yield estimation methods as processes occurring at each stage ultimately determine the final grain yield. There are typically three stages of rice growth with agronomic significance, namely vegetative, reproductive, and ripening, each with sub-stages or phases ( Figure 2) (Yoshida 1981). The vegetative stage is marked by leaf development and elongation of the main stem and auxiliary stems called tillers. The reproductive stage is divided into panicle initiation, booting, and heading phases. The panicle, the rice's flowering cluster, develops at the beginning of reproduction and is concealed initially. Then, booting occurs when bulging at the leaf sheath is detected. Heading is when the panicle protrudes above most of the plant's tillers, after which fertilization is completed within five days. At the ripening stage, grain weight increases and changes texture and color (milky, dough, yellow-ripe, and maturity stages) (GRiSP 2013). A rice plant completes its lifecycle in 90 to 160 days, depending on cultivar and environmental conditions (De Datta 1981).
Biomass accumulation and final yield are dependent on cultivar, environment, and management during each stage in the rice lifecycle (Horie, Yajima, and Nakagawa 1992). Breeders developed thousands of cultivars with selective yield-enhancing traits in different climate conditions (Khush 1997). Climate factors, such as temperature, sunshine hours, and precipitation, are critical determinants of rice yield with variable importance to the rice lifecycle. Each growth stage Figure 1. Diagram of concepts discussed in this paper on mapping rice yield using remote sensing. has its optimal temperature ranges. Deviation from the optimal conditions at any stage reduces yield significantly, even resulting in crop failure (Yoshida 1981;Li et al. 2015). Shading or limited incident solar radiation can also lower the yield (Yoshida and Parao 1976). Either deficiency or overabundance of rainfall, associated with extreme climatic events, can also decrease the yield (GRiSP 2013).
Remote sensors easily capture various physical changes associated with rice growth (Figure 3). Paddy rice is a particular type of crop requiring abundant water throughout its lifecycle except at maturity (Dong and Xiao 2016). The cycle of flooding, canopy development, draining, and drying in a paddy field is a distinctive feature from other crops. Most stages are physically distinct from one another. Careful selection of the spectral bands is critical to differentiate between stages and for successful mapping of paddy fields. The booting phase has the most suitable canopy for optical sensing Chang, Shen, and Lo 2005). At earlier stages of growth, there is minimal canopy cover, so background soil and water overwhelm the spectral response of rice plants. As the plant grows, leaves and stems contribute more to the signal until a saturation point, typically at maximum photosynthetic activity during the booting phase. The yellowing panicle and senescent leaves and heading are the best times for sensing rice yield ).
Although the single-date imagery can produce satisfactory results, multi-date images, such as those spreading over different stages, can improve the estimation accuracy . For instance, radar backscatter values are low immediately after sowing when paddy fields are flooded; then, it increases significantly during the vegetative stage ( Figure 4). From ripening to harvesting, radar backscatter decreases as canopy Figure 2. Agronomic growth stages of a rice plant (Oryza sativa). Main growth stages are vegetative, reproductive, and ripening with specific phases within each stage. Depending on the cultivar, it takes 90-160 days for rice to complete its life cycle from germination (of direct-seeded rice) or transplantation up to harvesting (source: adapted from GRiSP (2013)).

Figure 3.
Sentinel-2 Normalized Difference Vegetation Index (NDVI) and false color time-series imagery for 2019. During rice canopy development, spectral changes are easily detected from Sentinel-2 bands. NDVI is calculated from bands 4 (red) and 8 (NIR), depicted here as an increasingly dark green color at peak photosynthetic activity during flowering stages. Then, it decreases after maturity and during harvesting. In the false color composite of bands 3 (green), 4 (red), and 8 (NIR), greater vegetative cover and photosynthetic activity resulted in a deep red color. During harvesting, it changes to brown-gray. Pre-transplanting flooding was detected in rice plots encircled green and yellow. Plots encircled blue were harvested earlier, evident by its distinctive color in comparison to the surrounding plots. (Contains modified Copernicus Sentinel data 2019 processed by Sentinel Hub; use permitted under the Creative Commons Attribution License CC BY 4.0.). elements wilt (Chen and Mcnairn 2006). Mapping accuracy can be improved by using one image at the early vegetative stage and another at the reproductive stage (Shao et al. 2001).
In summary, the spectral response of rice plants and their environment varies during their lifecycle, creating useful occasions for accurate areal mapping and yield estimation of rice by means of multitemporal sensing. When using the single-date imagery, the best time to capture data is the booting phase. For the multi-date sensing, images at the early vegetative stage and reproductive stage can achieve the best estimation accuracies.

Optical properties
Structural and biological characteristics of rice influence its spectral response. Plant pigments, such as chlorophyll a and b and carotenoids, strongly absorb blue (450-485 nm) and red (625-740 nm) light while strongly reflecting near-infrared (750-1400 nm) radiation. Green light (500-565 nm) is also reflected, although much less than near-infrared, giving photosynthetic portions of leaves their green appearance ( Figure 5) (Woolley 1971). Aside from canopy features, the spectral signature of flooded paddy is affected by background water, suspended sediments, and soil, particularly at the earlier stages of growth (Martin and Heilman 1986). Higher amounts of soil nitrogen also result in greater nearinfrared reflectance values (Nguyen and Lee 2006). Measurable differences between phenological stages are discernable from field hyperspectral data. As canopy cover expands at later phases of growth such as booting and heading, the background influence diminishes, and eventually disappears at maturity. Regardless of the growth phase, rice (like most plants) always has the highest reflectance at near-infrared wavelengths ( Figure 5) while the most substantial increase and decrease in reflectance occurs at the red edge (700-740 nm) (Gnyp et al. 2014).
In selecting spectral bands, those maximizing the spectral disparity between rice and other crops should be used. Aside from individual bands, multispectral bands may be used to derive Vegetation Indices (VIs) to depict rice conditions better. The most widely used Normalized Difference Vegetation Index (NDVI) is given by Equation (1) below: where NIR is the near-infrared band reflectance and R is the red band reflectance value. NDVI is good at measuring plant vigor or greenness with a value ranging from -1.0 to 1.0. A higher NDVI value indicates . This irrigated rice area shows how detected backscatter changes over the growing season. Radar backscatter is the reflected signal detected by its source. Low backscattering is observed before planting due to flooded fields. Backscatter eventually increases due to scattering within the rice canopy and reflections from the water surface, peaking around the heading phase. At ripening, water content and canopy density decreases and so does backscattering. In between planting seasons, soil tillage increases surface roughness and hence increase backscatter before commencing flooding for the next cropping season. The red outline shows a contiguous are of rice plots, bisected by an irrigation canal. (Contains modified Copernicus Sentinel data 2019 processed by Sentinel Hub; use permitted under the Creative Commons Attribution License CC BY 4.0.).

Figure 5.
Spectral response curves different growth phases of rice. Due to changing distribution of canopy-soil-water cover in a pixel as the plant develops and increases in photosynthetic activity, reflected spectra changes from vegetative (tillering and stem elongation) and reproductive (booting and heading) growth stages, noticeably at the near-infrared wavelengths (source: adapted from Gnyp et al. (2014)).
more biomass or a healthier plant (Tucker 1979). NDVI is strongly correlated with photosynthetically active biomass, chlorophyll abundance, and energy absorption (Myneni et al. 1995). NDVI also corrects for various external influences such as topographic shadow and sensor noise. However, it is still liable to radiometric degradation in red and NIR bands, atmospheric effects, and background influence. Alternative VIs such as the Enhanced Vegetation Index (EVI) and the Soil Adjusted Vegetation Index (SAVI) can reduce signal degradation by using several empirical coefficients or additional bands such as the blue band (B) for EVI (Huete 1988;Huete et al. 2002). Their equations are: Substitution of the red band with red-edge in NDVI has considerably improved crop yield estimates . These alternative VIs can be derived from field hyperspectral data; thus, field spectral data provide a useful point of comparison with satellite image-derived results before scaling up ).

Radar properties
Persistent cloud cover during the growing season in major rice-producing regions such as Southeast Asia limits the usefulness of optical data. This led to the reliance on radar sensing, particularly Synthetic Aperture Radar (SAR). As an active system, radar has the "all-weather" capability, and can easily penetrate clouds and haze, thus providing continuous observation of paddies at all critical stages of growth (Inoue et al. 2002). The interaction of microwave signals with the target surface is measured as the backscatter coefficient (σ °), the signal reflected back to its source. It is sensitive to vegetation canopy structure and water presence of the underlying ground surface (Bouman 1995;Steele-Dunne et al. 2017). In rice paddies, total scattering is composed of volume scattering of canopies, multiple scattering between canopy elements and ground surface, and surface scattering from underlying soil and water . Studies have revealed a large dynamic range of σ° from the sowing to the maturity stages, and a strong relationship of σ° with rice biomass and height (Le Toan et al. 1997).
The choice of frequency, polarization, and incidence angles significantly influences σ°. σ° from lower frequencies (i.e. C-to L-bands) is strongly correlated with rice LAI and aboveground biomass due to deeper penetration of its canopy (Inoue et al. 2002;Clauss et al. 2018;Jia et al. 2014). Higher frequencies (i.e. X-, Ka-, Ku-bands) are more sensitive to changes in rice LAI at earlier stages (e.g. new transplantation), but the sensitivity is lost thereafter. C-, S-, and L-bands remained sensitive during the most life cycle of rice until the flowering stage (Inoue et al. 2002;Jia et al. 2014).
Radar polarization also influences σ°. HH (horizontal transmission and reception) generally has higher σ° values than VV (vertical transmission and reception), because of the attenuation of radar waves by vertical elements of the canopy, such as stems and leaves (Kuenzer and Knauer 2013;Tan et al. 2015;Bouvet, Le Toan, and Lam-Dao 2009;Nelson et al. 2014). However, HV (horizontal transmission, vertical reception) and VH (vertical transmission and horizontal reception) have equivalent σ. Modern radar sensors have fullpolarization modes, vastly improving the information available for monitoring paddies over older generations of sensors that are singly polarized at HH or VV only.
Finally, the temporal changes of σ° are more dynamic at a larger incidence angle (i.e. 45° or 55°) with lower frequency bands. Rice biomass is strongly correlated with L-band HH at 45° (r 2 = 0.99) (Tan et al. 2015;Inoue et al. 2002). At lower incidence angles (25-35°), biomass values hardly vary with the growth stage. Radar backscattering is also influenced by the band-polarization combination. For instance, L-VV is more loosely correlated with biomass than L-HH (Inoue et al. 2002;Jia et al. 2014). Therefore, the ideal configuration for rice yield estimation using radar sensing is low frequency (i.e. C and L) bands, HH-or multi-polarization, and a large incidence angle.
As the rice plant grows, backscatter intensity varies from one stage to another, subject to the aforementioned factors. Prior to planting, paddy fields have a low σ° value due to specular reflection, and hence appear as dark areas in radar images (Chen and Mcnairn 2006). At this stage, it is easy to differentiate flooded and non-flooded fields as the latter tends to have a higher σ° and appears as gray. During the vegetative stage, canopy elements (i.e. tillers, leaves, stems) have a higher σ° due to the expanded canopy volume and density, appearing as gray. For highfrequency SAR, such as X-band, σ° reaches the peak at early-to mid-tillering stages. Correspondingly, for low-frequency SAR such as C-band, the peak occurs during heading when rice appears as bright on imagery (Inoue et al. 2002;Jia et al. 2014). During ripening, the plant's scattering intensities decrease steadily until harvest. These temporal changes in rice backscattering are significantly different from those of other land covers and can effectively discriminate paddies from non-paddies (Chen and Mcnairn 2006;Le Toan et al. 1997;Ribbes and Le Toan 1999).
Similar to spectral indices for optical sensors, several radar-based indices have proven effective to detect rice biophysical properties. As most sensors are now available in multiple polarization modes, simple ratios between two differing polarization modes can highlight the difference between rice and non-rice crops, particularly when used for temporal analyses. Bouvet, Le Toan, and Lam-Dao (2009) employed HH/VV for C-band ENVISAT/ASAR for mapping rice areas in the Mekong Delta while Yang, Shen, and Zhao (2012) used the same ratio with C-band RADARSAT-2 for mapping rice extent using the ratio to assimilate parameters into a crop model. Frequency ratios of C-band VV and L-band VV from multiple-frequency data, such as those from field scatterometers, are negatively correlated with plant age as demonstrated by Oh et al. (2009). Data from the quad-polarization sensors allow the derivation of the radar vegetation index (RVI) which bears a strong correlation with leaf area index of rice (Kim et al. 2012;Kim and Van Zyl 2009). Its equation is given as: where σ vh is cross-polarized σ° values, while σ hh and σ vv are co-polarized σ° values. However, quadpolarization satellite data is scarce so another study a dual-polarized index, Dual Polarization SAR Vegetation Index (DPSVI), with Sentinel-1 data which was found to be strongly correlated with NDVI and aboveground biomass (Periasamy 2018). Its equation is given as: where σ vh is the cross-polarized σ° values of the SAR image, σ vv is the co-polarized σ° values, and max(σ vv ) is the maximum co-polarized σ° value in the image. These recent developments in SAR-based vegetation indices will be valuable for monitoring rice in cloudy regions where optical data cannot function properly.

Image resolutions
Rice cultivation area is an integral component of any rice yield estimation. Rice cultivation area derived from remotely sensed imagery is affected by its spatial resolution or pixel size. It significantly affects paddies of a small and patchy size, as in tropical Asia (e.g. <1 ha) (Lowder, Skoet, and Singh 2014). RS platforms are available in various spatial resolutions and can be chosen depending on the user's objectives. Use of high-or very high-resolution images at less than 1 m can minimize the mixed pixels of rice and non-rice. Unmanned aerial vehicles (UAVs)derived images provide fine spatial resolutions but tend to have a limited spectral range. Commercial satellites such as WorldView-3 and SkySat also feature sub-meter spatial resolutions with additional infrared bands. Low-resolution satellites (>250 m) such as MODIS and NOAA-AVHRR were useful for yield estimation over extensive areas, primarily highly uniform rice planting areas (i.e. irrigated paddies), but not for small paddies Son et al. 2014;Setiyono et al. 2018).
Recently, the proliferation of medium-resolution platforms (10 to 30 m) such as Sentinel and Landsat have paved the way for better detection of smaller farms (Siyal, Dempewolf, and Becker-Reshef 2015;Gilardelli et al. 2019; Multi-temporal images can vastly improve area and yield estimation, compensating for low spatial resolutions from other sensors (Basso and Liu 2019;Prathumchai et al. 2018). The aforementioned temporal developments of rice plant in section 2.1 can be utilized to effectively differentiate rice from other land covers such as urban, barren, or water bodies (Atzberger and Eilers 2011;Chen et al. 2011). These are better than the single-date images only when they are captured at different critical growth stages mentioned in section 2.1. A convenient method of acquiring multi-temporal images of the same paddy is to make use of drones or airplanes by timing image acquisitions to coincide with critical rice stages. Lowresolution platforms have the advantage of large swaths and revisit times short enough to observe critical stages of rice. However, these optical images can be contaminated by clouds. Additional problems may arise with optical imagery when it is necessary to correct for aberrant acquisitions of multi-dated images using temporal interpolation techniques (Chen et al. 2011;De Castro et al. 2018;Boschetti et al. 2009). These limitations can be addressed with the joint use of satellite SAR with optical images (Zhang et al. 2018). Table 1 presents several studies using various optical and radar platforms for yield estimation.

Key biophysical variables
Biophysical variables that reflect plant condition are reliable predictors of yield (Basso and Liu 2019). Retrieving these variables from RS imagery is a crucial prerequisite before being able to estimate yield. Both optical and radar imagery have been used to derive these variables such as LAI, biomass, Fraction of Absorbed Photosynthetic Radiation (FAPR), and Canopy Cover (CC) of rice. Crop models use these as primary inputs to estimate yield. This section discusses various approaches for biophysical variable retrieval, most of which are used as the inputs or canopy state variables in yield estimation models.
LAI is a dimensionless parameter of the one-sided green leaf area per unit of vegetated ground cover. LAI represents the available surface area that plants have to exchange components with the environment, including light, water and CO 2 (Carlson and Ripley 1997;Jonckheere et al. 2004). LAI value ranges from 0 to 10 with a higher value indicating a larger leaf area per unit ground area. It is a crucial variable in many crop models Keating et al. 2003). LAI is retrieved from VIs using either empirical relationships with ground values, physical model inversion, or nonparametric methods (Wang et al. 2016;Liang et al. 2015). Machine learning models (i.e. neural networks and support vector machines) estimate LAI from VIs with good accuracy if the image is acquired before the heading phase (R 2 = 0.9) (Jing Wang et al. 2016). In general, though the LAI-VI relationship is non-linear and crop-specific (Kang et al. 2016;Carlson and Ripley 1997). For radar, several scatterometer campaigns showed strong empirical relationships of σ° with LAI, particularly with the lower frequency L-, S-, and C-bands at a relatively large incidence angle with the HH-polarization (Tan et al. 2015;Jia et al. 2014;Inoue et al. 2002). σ° generally increases linearly with sensor configuration, and reaches a saturation point, typically during heading, but peaks can occur earlier with X-band (Inoue et al. 2002).
Biomass, the dry weight of rice organic material, is another variable studied from RS data. For optical systems, it is represented as a function of an accumulative metric over the growing period, usually cumulative NDVI (Tucker et al. 1985). RS can only Detect Aboveground Biomass (AGB), consisting of the total dry weight of foliage, shoots, reproductive organs, and grains. Other measurements related to biomass are Gross Primary Production (GPP) -the rate of organic material production by plants -and Net Primary Production (NPP) -the rate of GPP minus the rate of energy used for metabolic processes through respiration (Roxburgh et al. 2005). GPP and NPP are measures of rates and have a temporal dimension, while biomass refers to the actual weight of plant components at a specific time. FPAR, a measure of photosynthetic activity of the plant, has also been used to determine biomass accumulation in rice over time (Casanova, Epema, and Goudriaan 1998). FPAR, biomass, NPP, and GPP are strongly related to vegetation indices as demonstrated by empirical studies on agricultural areas (Tucker and Sellers 1986;Myneni and Williams 1994). The strong relationships with radar σ° were also established for biomass and plant height. Biomass is most strongly correlated with L, S, and C bands, but not at high frequencies (Inoue et al. 2002;Jia et al. 2014). The correlations are strong for plant height with HH-polarization, but not with VV polarization (Tan et al. 2015). Furthermore, C-VH σ° is highly correlated to FPAR, thus can derive measurements comparable to optical NDVI (Inoue, Sakaiya, and Wang 2014).
Canopy cover (CC) is another useful state variable with strong empirical relationships between RS data and field measurements. Also known as Fractional Vegetation Cover (FVC), it is related to LAI, both measuring canopy status. However, it addresses the saturation problem of LAI by representing green leaf area as a percentage cover over the ground area, as observed at nadir (Carlson and Ripley 1997). For rice agriculture, it is useful for monitoring light interception and crop water use, thus has been used successfully for yield estimation purposes (Prathumchai et al. 2018). Horticultural studies also showed good linear relationships (R 2 = 0.95) between CC and VIs, but they are untested for rice (Trout, Johnson, and Gartung 2008).
Aside from biophysical parameters, phenological parameters, such as the Start Of Season (SOS) and the Peak Of Season (POS) dates, served as reliable inputs into rice models to improve yield estimates. Deriving rice transplanting dates from MODIS time series data improved yield modeling outputs . This approach lessens the computational requirements as fewer images are needed in comparison to directly assimilating biophysical variables into the model.

Mapping rice cultivation areas
Areal measurements of rice fields are important to total rice yield that depends directly upon the area of cultivation. Rice cultivation area is commonly mapped from images based on automatic classification of pixel values. Several review papers have covered this topic comprehensively (Dong and Xiao 2016;Kuenzer and Knauer 2013;Mosleh, Hassan, and Chowdhury 2015). Table 2 provides a summary of traditional and new classification methods along with their advantages and disadvantages. Here, we discuss several new developments in mapping rice area. Pixel-based classification is the dominant and traditional method for image analysis; however, this is typically plagued by "salt-and-pepper" artifacts from poor algorithm performance, among other limitations. To address this, an approach called Geographic Object-Based Image Analysis (GEOBIA) groups pixels into clusters or segments of objects according to defined similarity criteria is increasing in popularity (Blaschke 2010;Blaschke et al. 2014;Hay and Castilla 2008;Kucharczyk et al. 2020). Rice paddies are highly fragmented spatially in some regions, but many of their physical features are homogenous. GEOBIA uses metrics such as shape and texture. Shape is the outline (geometry) of an object such as raised soil section of paddies while texture is the spatial variation of pixel values within a specific window (Gao 2009;Dong and Xiao 2016). These geometrical and topological features, in addition to spectral information, are key to the appeal of GEOBIA over per-pixel analysis. Recent studies using medium-resolution data with GEOBIA produced good results, for example, a study using both Landsat and RADARSAT-2 data achieved 86% accuracy by optimizing segmentation parameters and using a rule-based classifier to distinguish between rice paddies and permanently flooded water bodies (Pradhan, Sameen, and Kalantar 2017). Since GEOBIA can use various data sources, some studies extracted phenological information or used secondary vector datasets as inputs to improve segmentation results. These approaches obtained accuracies higher than 90%, outperforming pixel-based methods (De Castro et al. 2018;Zhang and Lin 2019;Singha and Sarmah 2019). Despite these advantages, compared to per-pixel methods, GEOBIA can be computationally costly and users assign segmentation parameters and, at times, decision rules for classification; hence, the optimal parameters must be tested for accurate results.
Sensor fusion, the synergistic use of multiple sensors, is used to combine the strengths of sensors in terms of spatial resolution, revisit times, or cloud penetration capabilities (Zhang 2010;Van Tricht et al. 2018). For example, STARFM synthesizes a dataset with Landsat's 30-m resolution at MODIS's higher revisit times (Gao et al. 2015). A study using this dataset showed a 6% improvement over Landsat-only methods . New efforts focus now on complimenting or fusing optical and SAR data. Various combination schemes of time-series Landsat, PALSAR, and RADARSAT data yielded accuracies up to 98% . Sentinel-1 and Sentinel-2 have also been used in complimentary ways. Using water and vegetation masks derived from Sentinel-2 after classification of multi-temporal Sentinel-1 resulted into 35% increase into non-paddy area accuracy (Inoue, Ito, and Yonezawa 2020). A scheme using random forests to combine both resulted into a 6% accuracy increase from only optical images (Steinhausen et al. 2018). These studies typically employ several schemes to combine satellite data hence careful testing of which imagery, bands, data products or derivatives to use will provide optimal results.
Recent improvements in classification algorithms focused on developing novel Machine Learning (ML) methods, primarily Deep Learning (DL). DL encompasses a wide range of technologies and applications and can be used in almost any step of remote sensing image processing but this topic is well beyond this paper so we focus on new developments for rice area delineation (see Ma et al. (2019) for an exhaustive review). DL, based on artificial neural networks, has demonstrated better performance than classical ML in combining spatial and temporal patterns in time-series data. Convolutional Neural Networks (CNN), Recurrent Neural Networks (RNN) and their variants are the most common types used. RNN in particular carries an appeal in phenology-based classification tasks due to its specialization for sequential data. Rice detection was improved  ML using CNN and RNN in France, Australia, Brazil, and China (Ndikumana et al. 2018;Guo, Jia, and Paull 2019;Filho et al. 2020;Zhao et al. 2019). In constructing these models, the results are dependent on optimal tuning of the parameters and architecture of the neural network along with a large amount of field reference data.

Models for rice yield estimation
Quantitative models are approximations of reality. Incorporation of observational data in models can improve simulation results considerably. Several types of models have been used to estimate rice yields spatially, including empirical, semi-empirical, and processbased crop models. Empirical models, also known as correlative or statistical models, are typically used over larger spatial scales such as the country or regional scale. They make use of long-term yield records and climatic trends (Lobell and Burke 2010). Process-based crop models, variously known as crop growth simulation models or simply crop models in the literature, are mathematical representations of the effects of weather, soil, water, stress, and management on plant growth and yield (Batchelor, Basso, and Paz 2002). Examples of multicrop models are APSIM (Keating et al. 2003), AquaCrop (Steduto et al. 2009), and DSSAT (Jones et al. 2003). Among rice-specific models are ORYZA , RiceGrow (Tang et al. 2009), and WARM (Confalonieri, Rosenmund, and Baruth 2009). Efforts to simplify components of process-oriented models produced semiempirical models, which use equations to partition simulated biomass into organ development over daily time-steps or represent scattering mechanisms occurring in the crop environment, resulting in total biomass and harvestable yield (Basso, Cammarano, and Carfagna 2013). Some purpose-built, remote sensing-based models were developed mostly from semiempirical models with RS variables (i.e. NDVI) as essential inputs .

Empirical models
Empirical models have been employed widely for trend analysis and to determine the relationships between several variables (Basso, Cammarano, and Carfagna 2013) due to increasingly data-rich resources from ground and satellite observations. They are relatively simple to use and easily scalable to homogeneous locations with the same cultivar. In-depth physiological knowledge is not required to construct and use these models. Through regression analysis, meteorological and even economic factors are used as inputs to these models. However, the established models are specific to areas where ground information is collected. Haphazard application of these models elsewhere seldom produces reliable results. It is rare for these models to be readily transferrable; hence, area-specific models are needed (Moulin, Bondeau, and Delécolle 1998). Empirical models have been used to estimate rice yields at a range of sites around the world (Table 3). A good example is hyperspectral data collected in the field. For instance, canopy state variables such as LAI under different nutrient, water, and cultivar conditions can be modeled from spectral data throughout the growing season (Maki and Homma 2014;Krishna et al. 2019). This modeling approach can establish if models are robust under different management conditions before being applied at the regional scale. Under ideal observation conditions, a good fit is achievable between the two, but in some cases, this approach can perform poorly (Wittamperuma et al. 2012). Rice yield can also be modeled at the field level empirically from different combinations of spectral indices derived from handheld spectrometers (Kanke et al. 2016). Temporal integration with other models can improve the reliability of underperforming models (Maki and Homma 2014).
Regression models of yield-radar backscatter show some correlation but these may be indirect as backscatter is more sensitive to parameters such as biomass, height, and geometry than to grain yield (Aschbacher et al. 1995). Furthermore, reliable radar-based yield estimates require at least three images -one from early-mid vegetative, at heading, and prior to harvest (Shao et al. 2001;Li et al. 2003). However, recent studies with crop cuts suggested that even at different polarization modes, linear relationships with backscatter are insignificant and other data sources or methods are more suitable (Guan et al. 2018).
Machine learning, a class of empirical modeling, provided new approaches to reliably estimate rice yields along with other types of crops (Palanivel and Surianarayanan 2019). Neural networks used with radar backscatter from RADARSAT-1 and yield measurements as inputs achieved an accuracy of 94% (Chen and Mcnairn 2006). Random forests were used with Sentinel-1 data with strong correlation (r 2 = 0.86-0.93) for three rice growing seasons of a single year with different sowing and harvesting dates (Clauss et al. 2018). Prediction of climate change impacts on rice yield is another application of machine learning.  used an artificial neural network with inputs of NDVI, temperature, solar radiation, and precipitation for yield estimation. Apart from better predictive performance (relative RMSE = 8-29%) than linear regression (13-60%), this approach is also flexible as inputs can be modified to include simulated climate variables under climate change scenarios as the authors have done. However, it involves meticulous selection of datasets and pre-processing steps. Like other empirical models, these may be useful at the regional or low-resolution level only (Gandhi, Petkar, and Armstrong 2016). This may result to less accurate results at smaller, specific scales if ground measurements are sparse or unavailable.

Process-based crop models
Process-based crop models tend to be more complex than empirical models, but they can better explain mechanisms of rice growth and, subsequently, yield. Due to the flexibility and relative simplicity of multicrop models, they have been commonly used by nonacademics in practical applications (Basso, Cammarano, and Carfagna 2013). Crop models can simulate the rice growth processes with considerable detail. Because of the detail of inputs required, they are suited only at point scales and are ill-suited for broad national-and regional-scale applications. Increasingly, RS data are used as direct inputs into crop models or to parametrize model runs to better reflect local conditions (Jin et al. 2018).

Rice crop models
With crop modeling technology maturing in the past couple of decades, various crop models have been developed. Those RS-based models are ORYZA (6 studies) and WOFOST (5), followed by SIMRIW (3), WARM (2), and DSSAT (2) ( Table 4). Although these models were not originally conceived for use with RS data, new models such as BESS-Rice and RS-P-YEC were built with imagery data or its derivatives as inputs.
ORYZA was developed by the International Rice Research Institute (IRRI) to model rice growth in optimal, nitrogen-limited, and water-limited conditions. Dry matter accumulation for plant organs is calculated daily as a function of the rate of production and the rate of phenological development. These are used to calculate dry matter at the end of the growing season (Bouman and van Laar 2006). This model was utilized with both optical and SAR data with LAI as the state variable, usually to estimate yield at a high accuracy (Huang et al. 2002;Shen et al. 2009;Guo et al. 2012;Yang, Shen, and Zhao 2012;Pazhanivelan et al. 2015;Setiyono et al. 2018).
The World Food Study (WOFOST) model was developed to simulate the growth of different crops in the tropics. It has become an integral part of Europe's crop yield forecasting framework owing to the general applicability of its biophysical core. Its main steps of simulation are phenological development, assimilation of carbon dioxide, transpiration, respiration, assimilated partitioning to different organs, and dry matter formation (De Wit et al. 2019). WOFOST was primarily used with satellite and field optical spectrometers to estimate yield (Kamthonkiat et al. 2010;Wu et al. 2013;Liu et al. 2014;Jin et al. 2015;Zhou, Liu, and Liu 2019).
Simulation Model for Rice-Weather Relations (SIMRIW) is a simplified process model that simulates the potential growth of irrigated paddies in relation to temperature, solar radiation and atmospheric CO 2 (Horie 1987). It is specific to rice and has been used for  (Park, Das, and Park 2018) Taiwan, China Random forest, support vector machine Heading date Yield MODIS NDVI Regional (Son et al. 2020) Abbreviations: LAI -leaf area index, NDVI -normalized difference vegetation index, EVI -enhanced vegetation index, SAVI -soil-adjusted vegetation index, OSAVI -optimized soil-adjusted vegetation index, MTVI2 -modified triangular vegetation index, WDRVI -wide dynamic range vegetation index, SR -simple ratio, TIPS -time-series index of plant structure, AGB -aboveground biomass, CC -canopy cover, LSWI -land surface water index.  Abbreviations: DM -dry matter, LAI -leaf area index, DVI -developmental index, RUE -radiation use efficiency, LNA -leaf nitrogen accumulation, SOS -start of season, RLGR -relative leaf growth rate, CC -canopy cover, IDVS -initial phenology development stage, TD -transplanting date, DWT -dry weight, ETo -evapotranspiration, SM -soil moisture, WRT -dry weight of roots, FPAR -fraction of absorbed photosynthetically active radiation, NDVI -normalized difference vegetation index, SAVI -soil-adjusted vegetation index, EVI -enhanced vegetation index, GNDVI -green NDVI, SCE-UA -shuffled complex evolution, PSO -particle swarm optimization climate change studies in rice-producing countries in Asia (Matthews et al. 1995). Its simulation procedure consists of modeling crop development stages, dry matter production, yield formation, and climatic water balance. Several studies used optical indices or derivatives to estimate yield or dry weight (Inoue, Moran, and Horie 1998;He et al. 2010;Raksapatcharawong et al. 2020). It has also been used with COSMO-SkyMed backscatter data to predict LAI with promising results after calibration (Maki et al. 2017). Developed in Italy, WARM improves upon previous crop models by introducing modules on peculiarities of paddy fields, such as micrometeorology, soil hydrology, cold shocks, and blast diseases. Additional modules on soil hydrology improved on temperature effects of floodwater during irrigation (Confalonieri, Rosenmund, and Baruth 2009). Studies with optical data resulted in highly accurate estimates at the farm and regional levels in rice growing areas in Europe, with LAI being the state variable used to calibrate the model (Gilardelli et al. 2019;Pagani et al. 2019).
Decision Support System for Agrotechnology Transfer (DSSAT) has various modules for different crops. Its core comprises a collection of multiple independent crop models that describe crop parameters with similar soil modules, aiming to have a standardized framework for simulation as it was designed for wide adoption beyond research (Jones et al. 2003). For rice, the core module is CERES-Rice, which considers inputs and effects of weather, genetics, soil water, soil carbon and nitrogen content, and management factors. Single or multiple seasons and crop rotations are also considered (Singh, Singh, and Hasan 2007). The model with optical satellite data significantly improved the estimated yield at the regional scale (Son et al. 2016;Rezaei et al. 2016) BESS-Rice is the rice variant of the BESS, which uses an enzyme kinetic model to estimate GPP, a more mechanistic method than empirical approaches, such as LUE-based models . BESS-Rice improves upon this by introducing a process-based dry matter allocation module. It has reliable performance for rice yield estimation, with a RMSE of 534.8 kg· ha -1 (7.7%) and a bias of 242.1 kg· ha −1 (3.5%).
The RS-P-YEC model is based on the BEPS (the Boreal Ecosystem Productivity Simulator) model which simulates terrestrial NPP. BEPS was improved by adding a leaf canopy physical module that modeled energy conversion into biomass. RS-P-YEC was initially used for winter wheat yield (Wang, Yang, and Yagola 2011). Fed with MODIS LAI data, this model produced a relative error of 7-9% in estimating rice yield in several sites across six provinces in China ).

Data assimilation
Data assimilation aims to incorporate canopy state variables such as LAI and CC with RS variables such as NDVI to optimize model parameters (Jin et al. 2018;Huang et al. 2019). These different types of data can be assimilated via three methods of calibration, forcing, and updating. RS data have been assimilated with crop models using different combinations of platforms, crop models, and algorithms. Table 4 summarizes common state variables considered, RS inputs, and model outputs at various locations and scales around the world.
Calibration entails the use of an optimization algorithm to adjust initial parameters of the model to obtain optimal consistency between simulated values and imagery-observed values. The crop model is run manually or automatically with model parameters set within realistic ranges. Commonly used algorithms to minimize differences between observed and modeled values are simple search algorithms (Maki et al. 2017;Pagani et al. 2019;Gilardelli et al. 2019), shuffled complex evolution Yang, Shen, and Zhao 2012;Shen et al. 2009), and particle swarm optimization (Son et al. 2016;Wang et al. 2014;G. Zhou, Liu, and Liu 2019). Algorithm selection depends mainly on computing power and compatibility with the crop model being used. For instance, a study assimilating rice phenology information and EVI into WOFOST resulted in a computation time of 3.8-13.7 h . Careful selection of RS products or derivatives can shorten computational time while improving the accuracy of results.
Forcing refers to the direct replacement of simulated values with RS variables that must have the same frequency as the temporal increment of the crop model. This requirement is quite challenging to fulfil with certain types of satellite data. Strictly speaking, forcing methods are not assimilation strategies as they only replace simulated canopy state variables. The introduction of NDVI measurements to the ORYZA1 model achieved satisfactory accuracies at the regional level (relative error: -6 -6%) (Huang et al. 2001). Similarly, replacement of CC measurements in the AquaCrop model by NDVI achieved a Mean Absolute Error (MAE), RMSE, and correlation coefficient (R 2 ) of 0.159, 0.182, and 0.882, respectively (Prathumchai et al. 2018).
Updating involves sequentially updating model state variables whenever new observations are made available. Underpinning this approach is the assumption that better simulation data at day t will improve the accuracy of the simulation data at succeeding days. This approach is more flexible and computationally faster than both forcing and calibrating methods. However, it is more computationally expensive with a high measurement uncertainty (Jin et al. 2018). This method of assimilation has not been used particularly for rice, despite its promising performance with other crops (Dorigo et al. 2007;Jin et al. 2018). As computational power is expected to improve in the future, along with advances in RS platforms, updating methods may be used much more widely in crop model data assimilation.

Semi-empirical models
Midway between purely statistical and crop models are semi-empirical models. They retain the discrete time steps of crop models but parameterize many of the equations using experimentally derived information (Delécolle et al. 1992). Monteith's efficiency concept is a prime example of semi-empirical models, where healthy plants are assumed to absorb radiation at a Light Use Efficiency (LUE) conversion factor in producing AGB or dry mass (Monteith 1972(Monteith , 1977. This factor is typically assumed to be a constant specific to vegetation or crop type under observation but can be determined experimentally (Mitchell, Sheehy, and Woodward 1998). The relationship of AGB to LUE and PAR is expressed as: where FPAR is the fraction of absorbed PAR (APAR) to total incoming PAR. This variable can be estimated from RS images based on its close relationships with VIs (Clevers 1997;Myneni and Williams 1994). PAR can be measured from meteorological stations or can be assumed to be approximately half of total incoming radiation. FPAR is also related to LAI using Beer-Lambert's law (Boschetti et al. 2011). LUE has been routinely used to provide regionwide rice yield estimates to varying degrees of success. Guan et al. (2018) fused Landsat-MODIS EVI product to estimate rice yield in Vietnam at a moderate accuracy (relative RMSE: 17-19%). In another study, it was used to predict NPP of rice in the Continental United States at a relatively poor accuracy (R 2 = 0.53, RMSE = 3.42 t/ha) (Marshall, Tu, and Brown 2018). The benefit of using semiempirical modeling is the ease of scaling the estimation to region-wide yield. However, the reliance on optical imagery reduces the estimation accuracy as specific data crucial to calculate LUE are missing. Another semi-empirical method is to dynamically model Harvest Index (HI), the fraction of harvestable grains to total AGB (Boschetti et al. 2011). In the modeling, HI is usually assumed to be constant for a given type of crop, but is affected by climatic and soil stress, both of which vary with time.
So far, a few rice-specific semi-empirical RS-based crop models have been developed, including GRAMI (Table 5). GRAMI is one of the earliest RS-based semi-empirical models developed by Maas (1993aMaas ( , 1993bMaas ( , 1992. A modified version has been used for rice crop in South Korea using UAV data (Kim et al. 2017;. It also has a within-season updating module, which can be considered as a type of updating data assimilation tool . Common to these models is that RS variables are used either as a canopy state variable directly or to update model state variables when observations are available. In essence, all these semiempirical models encompass the LUE concept as its core. Typically, they do not have water balance and soil balance, so semi-empirical models cannot show how changing environmental conditions influence rice yield. If estimated from radar data, yield is retrieved by inverting semi-empirical scattering models (Table 5). These models simulate the complex scattering mechanisms occurring within canopy elements and between rice and soil surface. At the "forward model", their inputs consist of ground measurements and estimate backscatter values ). These were found to be useful to interpret field measurements but were not fully reliable for variable retrieval due to their inability to accurately simulate scattering within a canopy, the large number of required inputs, and inherent limitations within the model itself (Inoue, Sakaiya, and Wang 2014). In general, inversion continues to be a significant challenge in radar-based retrievals of biophysical variables and is a delicate process varying with sensor configuration and time of observation, but improving the optimization method can potentially produce reliable retrievals (Wang et al. 2011;Zhang et al. 2017). Promising results with the Genetic Algorithm (GA) optimization with the Rice Canopy Scattering Model (RCSM) produced low errors for two parameters (rice ear length and ear number density) which were used to map rice yield with some overestimation when compared to statistics . Despite several limitations of radar-based semi-empirical models, a potential approach to improve yield estimate is by integrating these models with process-based models (Inoue, Sakaiya, and Wang 2014).

Linking models to spatial information
The yield estimation models discussed earlier are upscaled from point-scale to regional-scale by executing the model at the pixel-level, resulting to high-granularity yield maps. Other gridded data such as climate data or other remote sensing derivatives are also vital elements as they serve as inputs to the model. Spatial representation of additional products such as LAI, biomass, plant age, and phenological parameters are also possible, which have value in assessing farm conditions over the growing period Campos-Taberner et al. 2017;Pagani et al. 2019). Government agencies typically collect yield statistics at administrative levels. Accuracies of yield estimation models are compared against these by aggregating pixel-level estimates into administrative-level averages, using mostly RMSE, MAE, and R 2 measurements as error or agreement measures (Gilardelli et al. 2019;Guan et al. 2018). Although these models can be readily linked to spatial data, there are several unresolved issues related to models themselves and other parameters in the models, such as uncertainties in paddy cultivation area and modeling approach, that can degrade the reliability of estimated rice yields.

Knowledge gaps
The use of RS with quantitative modeling techniques for rice extent mapping and yield estimation have numerous areas for improvements. Further improvements are possible owing to the maturity of crop modeling as a discipline which has since moved from a purely academic exercise to be an integral part of decision support systems and models at various spatial scales. Parallel to these developments is rapid advances in RS technology, from more satellite systems providing a great variety of products at higher frequencies to mainstream use of drones for precision agriculture. This section identifies critical knowledge gaps and areas of improvement in rice yield estimation.

Knowledge-and phenology-based mapping of paddy type
The image classification methods mentioned in Section 2.4, either pixel-based or GEOBIA, are good at classifying all types of paddy fields as a single land cover. On the other hand, precise spatial estimation of rice yield requires differentiation between irrigated, rain-fed, and upland paddy fields. Current methods relying on pixel properties are incapable of accomplishing this task. Hence, non-imagery data are needed to generalize typical characteristics of these rice ecosystems and to synergize research efforts with operational needs. Spatial variability and agronomic practices in major rice-producing countries render it challenging to use the information solely from digital imagery. Employing a knowledge-based method with various inputs such as topographical information, crop calendars, climate data, irrigation maps, and other locally relevant inputs can potentially improve existing methods to classify paddies. For instance, Digital Elevation Models (DEMs) are useful to differentiate between upland and lowland areas. Irrigation network maps can also provide additional information helpful to distinguish between irrigated and rainfed   Abbreviations: LAI -leaf area index, NPP -net primary productivity, NDVI -normalized difference vegetation index, AGB -aboveground biomass, RDVI -renormalized difference vegetation index, OSAVI -optimized soil-adjusted vegetation index, MTVI -modified triangular vegetation index, EVI -enhanced vegetation index, GCVI, green chlorophyll vegetation index paddies. Crop calendars can also refine estimation methods by providing general but localized schedules of planting, useful when crop rotation is practiced.
Phenology-based classification is an increasingly common method. Its strength for agricultural applications is due to the highly dynamic nature of croplands. In the case of rice, phenological stages and changes between flooding, canopy development, grain-filling, and harvests are easily distinguishable from less dynamic covers such as forests, water, and built-up areas. Their shorter durations also distinguish them from many perennials. In spite of this, annual crops, mainly cereals like corn, can be difficult to distinguish in crop rotation areas due to similar phenology. Additionally, rice ecosystem types can be challenging to detect using solely time-series spectral data. Utilizing additional spatial information and knowledge-based approaches in addition to phenological indicators can potentially achieve this.

Alternative canopy variables and next-generation sensors
LAI that can be reliably and easily acquired from RS data has been used widely in crop models. However, the inconsistent relationship of LAI with most VIs and the saturation problem causes uncertainties to yield estimates. Although LAI is the primary canopy variable in crop models, other canopy variables, such as CC have also been used with the AquaCrop model (Steduto et al. 2009). The model, formulated by the UN-Food and Agriculture Organization is built on the concept of water productivity -the proportion of yield to evapotranspiration, indicating water use efficiency to produce biomass (Steduto, Hsiao, and Fereres 2007;Brady 1981). It makes use of CC instead of LAI (Steduto et al. 2009;Silvestro et al. 2017). Compared to LAI, CC is easier to derive from RS as it does not saturate and consolidates leaf structural attributes (Silvestro et al. 2017;Prathumchai et al. 2018;Steduto et al. 2009). CC can be derived from the empirical relationships between RS data and ground measurements using pixel decomposition, machine learning, or inverting physical models like that of LAI .
Satellite imagery has proved popular in monitoring rice growth in the past, owing to its relatively low cost, frequent revisits, and broad areal coverage suitable for region-wide monitoring. However, at smaller scales and in precision farming applications, space-borne sensors can be inadequate. As costs associated with drones and their payloads are decreasing, they are more useful in studies that require fine resolution images obtained at a high frequency. This review covered a few drone-based studies with optical payloads, but potentials with active microwave sensors such as radar and LiDAR have been demonstrated with non-rice crops, providing highprecision biomass estimates derived from canopy height (Moeckel et al. 2018). New spaceborne LiDAR sensors, such as NASA's ICESat-2 and Global Ecosystem Dynamics Investigation (GEDI) and JAXA's forthcoming Multi-footprint Observation Lidar and Imager (MOLI), will further expand the possibilities of LiDAR applications to vegetation monitoring (Salas 2020). Another recent application showed that combining deep learning algorithms with multi-source imagery including LiDAR and optical sensors significantly improves crop detection (Prins and Van Niekerk 2020). Drone data are also useful for studying breeding activities by providing high throughput phenotyping as an alternative to traditionally laborious, time-consuming phenotyping methods Tattaris, Reynolds, and Chapman 2016).

Modeling approaches and hybrid models
Some parameters in process-based crop models are treated as a constant despite their known dynamic nature. Among them, the most important is the harvest index, a determinant of final grain yield relative to AGB. It is also an integral component in semi-empirical models. Several crop models employ a constant, crop-specific HI which can disregard the negative influence of droughts and stressful temperatures at the reproductive stage on HI (Moriondo, Maselli, and Bindi 2007;Peng et al. 2004;Prasad et al. 2006). The use of a constant HI has resulted in less-than-ideal model performance (Guan et al. 2018;. Some of the discussed models, such as ORYZA and SIMRIW, do dynamically estimate HI and have been useful for climate impact assessment (Matthews et al. 1997). New models such as AquaCrop were specifically developed for water stress or irrigation applications directly influencing HI, and they present new opportunities to improve upon ORYZA and SIMRIW with potentials to be extended to non-rice crops as well.
Point-based model results can be upscaled to the entire study area or a region using several methods, one of which is the Scalable Crop Yield Mapper (SCYM) (Lobell et al. 2015). This model uses "pseudoobservations" by generating through Monte Carlo sensitivity simulation of several hundred runs of a crop model. LAI of the model can then be empirically related to a satellite-based index, such as the green chlorophyll VI, and then up-scaled to the regional extent using Landsat images. This modeling approach has been used for African corn and wheat farms, India's Wheat Belt and the US Corn Belt with reasonable accuracy (Azzari, Jain, and Lobell 2017;Jain et al. 2013;. However, the performance of the same novel modeling approach to different rice cultivars using ricespecific models has not been explored yet. Therefore, it remains unknown how much improvement can be expected in model accuracy. Recent advances in machine learning also present opportunities for hybrid modeling of rice yield. Researchers have improved predictions of wheat yields by using the APSIM model with two separate machine learning algorithms, Random Forest (RF) and Multiple Linear Regression (MLR) . The APSIM + RF model was able to explain 81% of the variation in the observed yields, a 33% improvement from APSIM alone, and 19% from APSIM + MLR. This novel method can be useful with remote sensing which can supply biomass and phenological data in the absence of field data. By taking advantage of the predictive power of machine learning and explanatory features of crop models, these modeling methods can improve yield predictions in data-rich areas.

Future rice yield
Global temperature is likely to warm by 1.5°C from the pre-industrial level in the next three decades likely to reduce crop yields globally (IPCC 2018). Past studies on climate change effects on rice yield vary widely in scope and scale but are significantly limited to the region-or nation-scale. Climate effects on crop yield have been observed with on-site experiments or simulated with crop models for point locations (Xu, Wu, and Ge 2018;Van Oort and Zwart 2018). Point-scale rice models, such as ORYZA and SIMRIW were used extensively in an early study involving multiple countries (Matthews et al. 1997(Matthews et al. , 1995. However, point-based crop models are logistically and practically challenging to upscale regionally due to extensive fieldwork and necessary model calibration. Region-wide studies with field-level specificity are increasingly needed as a decision support tool and for future food security assessments. Only  studied the impacts of climate change on national rice yields in South Korea, using artificial neural networks with a database of meteorological, rice yield, and satellite data. This approach is useful for data-rich areas with long-term records of climate and yield. It is less useful for regions with sporadic or patchy data. One solution is to combine long-term RS observations or their derivatives (i.e. phenological information) with crop models that can accept projected climate changes as inputs. Model sensitivity to temperature and rainfall changes can be tested to reveal potential yields under various warming scenarios. The simulated results may be upscaled to assess regional to global effects of climate change on food security and yield gaps (Lobell, Cassman, and Field 2009;van Ittersum et al. 2013) in the next 50 years.

Operationalization
As this research field continues to mature, operationalization of yield estimation is a logical way forward. An effective operational system lies on integrated approach that incorporates the strengths of both ground-based calibration activities and powerful computational infrastructure. Currently, such a system solely for rice exists only in the Philippines through the Philippine Rice Information System (PRISM, https://prism.philrice.gov. ph/) and a similar IRRI project, RIICE (https://www.riice. org/), implemented mostly in Vietnam, is currently ongoing. There are also potentials for incorporating yield estimation models into more applied and management-oriented products such as suitability mapping (Ujoh, Igbawua, and Ogidi Paul 2019), disease and yield loss assessment (Mfuka, Byamukama, and Zhang 2020), vulnerability analysis (Nistor 2019;Dela Torre et al. 2019), and mobile applications (Veerakachen and Raksapatcharawong 2020). Eventually, uptake of this technology is beneficial at all levels of stakeholders, from policymakers to regional agricultural managers, and ultimately, to rice farmers.

Conclusions
The parallel development of RS technology and quantitative methods to map rice extent and model crop growth variables has created a huge potential for estimating rice yields with an enhanced accuracy. This review comprehensively summarizes and critically assesses prevailing methods of estimating rice yields using quantitative modeling and geospatial tools. At present, only VIs, radar backscatter, LAI, NPP and FPAR are commonly derived from satellite data and used as inputs into crop models that can be empirical, semi-empirical, or process-based. More accurate modeling results may be attainable by using more novel variables, such as canopy cover or canopy height from LiDAR data. Crop models (WOFOST, ORYZA, GRAMI, and SIMRIW) were coupled with RS data through three methods of data assimilation (calibration, forcing, and updating). More purpose-specific models can be created by adding new modules sensitive to pests, diseases, or extreme events to existing models. However, only a limited number of models have been developed for rice. So precise yield estimation by rice cultivar has not been tried, even though it is possible to map rice into broad types (upland, irrigated, and rain-fed) using knowledge-based image classification. The current crude estimation (e.g. HI treated as a constant) can be made dynamic by assimilating climate information with rice crop models fed with RS derived crop parameters. With the use of these improved modeling approaches, it is expected that our current knowledge gap in regional rice yield in light of the projected climate change will be bridged soon.

Acknowledgments
This study was supported by the New Zealand Ministry of Foreign Affairs and Trade PhD Scholarship to the first author. The authors would also like to thank the Copernicus Programme and Sentinel Hub for Sentinel imagery. Lastly, we are grateful to the valuable and constructive comments of the two anonymous referees that have greatly