In this study, we considered the whole country of Bangladesh as our study area (Fig. 1). It lies between 20° 34′ to 26° 38′ North latitude and 88° 01′ to 92° 41′ East longitude (Fig. 1b). The country has an area of approximately 143,998 km² with a total population of 160 million (https://www.cia.gov/library/publications/the-world-factbook/geos/bg.html). The study area experiences almost flat topography except for the Chittagong Hilly areas in the southeastern part, the Sylhet mountain range in the northeast, the Madhupur Tract in the central part, and the Barind Tract in the northwest. In terms of climate, it observes sub-tropical conditions with three dominant seasons, i.e., dry winter (December to February), pre-monsoon hot summer (March to May), and the rainy monsoon season (June to October) (Rashid, 1991). Temperature-wise, January is the coldest month with approximately 17°C and April is the hottest one with temperature in between 30-40°C. On an average per year, Bangladesh receives a minimum of 2000 mm of rainfall over most of the regions except the northwestern region of Rajshahi (~1600 mm), and coastal area of Chittagong in the southeast and northeastern part of Sylhet (~4000 mm) (http://www.weatheronline.co.uk/reports/climate/Bangladesh.htm). It is interesting to mention that approximately 80% of the rainfall happens during the monsoon season (i.e., June to October) and negligible amount in the dry season (i.e., December to February). Despite some variability in climatic conditions (in terms of both temperature and rainfall regimes in particular) throughout Bangladesh, the boro cultivation primarily depends on the ground water based irrigation schemes, which is usually being controlled by the farmers. Thus, we considered not to investigate the impact of natural climatic differences on the boro production across the country. Also, note that Bangladesh is a natural disaster (that include flood, cyclone, drought, river erosion, etc.) prone country, which plays critical role in the agricultural productivity including boro rice.

Bangladesh is an agricultural-dominant country, where approximately 62.8% of the total area is dedicated towards crop production. Among all crops, rice is the major contributor, i.e., cultivated in approximately 75% of the agricultural land (Karmokar & Shitan, 2012). In general, there are mainly three different types of rice are cultivated, i.e., aus (during April/May to July/August), aman (during April/May to November/December), and boro (during January/February to May/June). It is worthwhile to note the boro rice is the dominant contributor [i.e., approximately 55% of the total rice production during 2012 season (BBS, 2012)] among all kinds of rice production.

Here, we employed MODIS-derived two 16-day composite of NDVI images at 250 m resolution (i.e., MOD13Q1 version 005) during the initial (i.e., January 1 to January 16) and peak greenness (i.e., March 23/24 to April 6/7 depending on leap year) stages for the period 2007-2012 associated with the boro rice developmental stages. These images were corrected both geometrically and atmospherically; and quantified in the scope of NASA’s earth observing system program (Huete et al., 2002; Justice et al., 2002). Upon acquisition of these images, we performed a set of pre-processing prior to the development of the rice forecasting method, such as (i) re-projecting of the NDVI images into Transverse Mercator from the default Sinusoidal Gridded Projection system; (ii) mosaicing of the two tiles of NDVI images in order to have continuous images for the entire study area; and (iii) identifying cloud-contaminated pixels using the quality assessment information, and excluded from further analysis.

Apart from the MODIS data, we employed two ancillary datasets, such as (i) boro suitability map available from Bangladesh Agricultural Research Council (BARC) to aid the delineation of the area under boro rice cultivation, and (ii) ground-based information about boro rice area and its yield available from BBS to calibrate and validate the performance of results obtained from MODIS-derived boro rice estimates of both area and its yields. The boro suitability map was developed at 1:250,000 scale as a function of soil permeability, soil depth, soil moisture, nutrient regimes, soil salinity, soil pH, soil consistency, drainage, inundation depth, flood hazards, slope, length of the crop growing season, and temperature regimes (Hussain et al., 2012). This map was released in 2012 and consisted of five categories, such as very suitable, suitable, moderately suitable, marginally suitable, and not suitable for boro cultivation. In case of ground-based information, they were collected using regular sampling surveys and interviewing the farmers.

The methods had two major components, (i) delineating areas under rice cultivation before harvesting; and (ii) forecasting rice yield as a function of NDVI and assessing the accuracies.

Upon applying the thresholds determined in sub-section “Delineating areas under rice cultivation”, we applied them over validation dataset during the period 2010-2012 and the results are shown in Fig. 5. It revealed strong relations (i.e., R ² in the range 0.93 to 0.95 and RMSE in the range 30,519 to 37,451 ha, at 23 district-levels) between modeled and ground-based area estimates. Our results were also similar to other studies: (i) Panigrahy et al. (1997) obtained accuracies in the range 85.8 to 91.5% in the early estimation of rice areas in Wet Bengal, India; (ii) Patel et al. (2004) found approximately 95% accuracy in forecasting Kharif rice acreage in Orissa, India; (iii) Chen et al. (2011) observed a R ² of approximately 0.76 while compared the remote sensing-based estimates with Bureau of Agricultural Statistical data in Nueva Ecija, Philippines; (iv) Guo et al. (2014) observed an overall accuracy of 97.3% in determining rice areas during early season in Southern China; (v) Guan et al. (2016) used MODIS-derived time series images of NDVI over the entire region of Vietnam and observed high mapping accuracy (R ² = 0.81) while comparing with the statistical data; and (vi) Jin et al. (2016) used Landsat-derived vegetation indices, such as NDVI, EVI, and land surface water index for mapping the rice area over the Sanjiang Plain of northeast China, and found very good results (R ² = 0.85) while comparing with the census data.

In addition, we also compared the boro rice area estimated using only two images acquired during the initial/transplantation and peak greenness stages in the scope of this paper with those extracted in an earlier study (Mosleh & Hassan, 2014). In that study, we employed ten 16-day composite of NDVI images acquired over the entire growing season and estimated the boro rice area after the harvesting of the crop. It would be interesting to note that the outcomes of the proposed method outperformed the earlier one. For example, the overall RE’s were in the range ±10% using the proposed method in comparison to about ±30% using the method described in the earlier study (Mosleh & Hassan, 2014) (see Fig. 6 for details). Our findings implied that more images induced more error similar to other studies, e.g., (LeGrand, 2002); which were most likely related with data/image redundancies. In order to address this, one of the possible alternates would be the employment of ‘principle component analysis’ over all the ten images acquired during the entire growing season (e.g., Mosleh & Hassan, 2014); which was beyond the scope of this paper. Note that both methods would complement each other as understanding the area (thus the yield as described previously) approximately one to two months prior to harvesting would be critical in ensuring food security. However, only the earlier estimation would hold if any natural disaster would not take place.

Figure 7 shows the relationship between the district-specific average NDVI values at the peak greenness period (i.e., in between March 23/24 to April 6/7 depending on leap year) and ground-based rice yield estimates during 2007-2009. It demonstrates strong relationships (i.e., R ² in the range of 0.75-0.89), which are similar to other studies (Nusara et al., 2011; Noureldin et al., 2013).

The relationship observed in Fig. 7d (i.e., three years, 2007-2009, of combined data) was implemented in forecasting rice yield during 2010-2012; and is shown in Fig. 8. It reveals strong relationships (i.e., R ², RMSE, were in the range 0.76 to 0.86, 0.21 to 0.29 Mton/ha, respectively during the period 2010-2012) between ground-based and forecasted yields at 23 district-levels. In addition, we also observed that REs were in the range -5.45 to +6.65% when compared at 23 district-levels (Fig. 9). In fact, our findings were similar to other studies, such as: (i) Nuarsa et al. (2011) found R ² ≈ 0.93 over Bali Province, Indonesia; (ii) Rahman et al. (2009; 2012) observed reasonable relationships (i.e., R ² ≈ 0.56 for aus rice and R ² ≈ 0.89 for aman rice) over Bangladesh; (iii) Chang (2012) reported good agreements (i.e., R ² in the range 0.57 to 0.61) over Shi-ko, Taiwan; (iv) Huang et al. (2013) predicted the rice yield over five rice growing provinces of China and observed good results (i.e., R ² in the range 0.84 to 0.97, and overall RE of 5.82%); (v) Noureldin et al. (2013) demonstrated strong relationship (i.e., R ² in the range 0.90 to 0.95) over Kafr El-Sheikh Governorate, Egypt; and (vi) Kar & Kumar (2014) reported good results (i.e., R ² value of 0.75 during flowering stage) over Orissa, India.

Despite good agreements, it would be worthwhile to note that our forecasting would hold if the rice crop not be affected by natural disturbances (that include cyclone, insect outbreak, etc.). In addition, approximately 14 to 24% of disagreements between the ground-based and forecasted rice yield estimates could be attributed by other factors, such as (i) satellite images might be affected by atmospheric effects (e.g., cloud), which degrade the quality of the acquired data and thus the developed crop-yield model (Mkhabela et al., 2011); (ii) variation in climatic conditions at micro-level during the growing season could potentially impact the agreement level of rice yield (Son et al., 2013; Mosleh et al., 2015); and (iii) uncertainty associated with ground-based yield estimates due to insufficient observations could lead to poor rice yield assessment (Mosleh & Hassan, 2014).

Although our proposed method demonstrated its effectiveness in forecasting rice yield; we would require to consider the following issues in order to potentially enhance its capability, such as (i) the use of relatively high resolution satellite imagery would able to distinguish rice cultivated area more precisely in particular to heterogeneous agricultural practices; (ii) it would be necessary to revisit the calibration parameters, i.e., slope and intercept of the expression illustrated in Fig. 7d periodically; and (iii) it would be worthwhile to investigate other important weather variables, such as temperature, precipitation, and solar radiation in particular; as they were also used in the of rice yield forecasting framework (Savin & Isaev, 2010). In addition, we strongly recommend that the proposed modelling schema be thoroughly evaluated prior to implementing in other geographical locations.