Development of a Remote Sensing-Based “Boro” Rice Mapping System

Abstract

Rice is one of the staple foods across the world, thus information about its production is essential for ensuring food security. Here, our objective was to develop a method for mapping “boro” rice (i.e., cultivated during the months January to May) in a Bangladeshi context. In this paper, we used a Moderate Resolution Imaging Spectroradiometer (MODIS)-derived 16-day composite of normalized difference vegetation index (NDVI) at 250 m spatial resolution in conjunction with ancillary datasets (i.e., land use map, crop calendar, and ground-based rice production information) during the period 2007–2012. The proposed method consisted of three procedures: (i) ISODATA clustering and determining the boro rice signatures in temporal dimension using data from the period 2007–2009; (ii) formulating a mathematical model for extracting the boro rice areas using data from the period 2007–2009; and (iii) model calibration using data from the period 2007–2009 and its validation using data from the period 2010–2012. The implementation of the abovementioned procedures revealed reasonable agreements between the model (i.e., MODIS-based) and ground-based estimates of boro rice area at both country (i.e., percentage error in the range −0.83–1.42%) and district levels (i.e., r² in the range 0.69–0.89) during the period 2010–2012. Our proposed method demonstrated its effectiveness in mapping rice system at the regional/country scale.

1. Introduction

Rice is one of the most important crops that support about 40% of the caloric demands for the people across the world [1]. Significant amounts of rice (i.e., greater than 80% of the global production) are usually produced in Asia [2]. Among the Asian countries, Bangladesh (i.e., one of the most densely-populated countries in the world having ∼160 million people with a total area of 143,998 square kilometer [3]) is the sixth largest rice producer at the world-scale [4]. However, the sufficient local production of rice encounters several challenges, such as (i) the increasing demands on rice due to an overly growing population (i.e., from ∼80 million to ∼148 million during the period 1980–2010 [5]); (ii) loss of rice producing lands of the country due to both anthropogenic interventions (e.g., developing infrastructures, housing, and industries) and natural degradations (e.g., erosion of the major rivers); (iii) increased frequency and magnitude of the several climate-induced natural hazards (that include monsoon/summer flooding, cyclones in the coastal areas, drought in the dry season especially in the northern region, and saline intrusion in the south-west region of the country) that significantly impact/affect rice production; and (iv) conversion of rice-producing lands into other farming systems (e.g., shrimp farming in the south-west region in particular). The magnitude of these challenges is always dynamic, thus it would be worthwhile to develop a rice mapping system that would be able to estimate total production in order to ensure food security for the people of Bangladesh.

One of the most common and widely used methods for estimating rice cultivated areas is the use of agricultural statistical data acquired through field visits and interviewing the farmers. Despite its invaluable ability for understanding historical trends in rice production, this method is extremely tedious, time-consuming, inconsistent, too generalized, and unable to provide information in a timely manner as well as delineate detailed spatial distribution of areas under rice cultivation. In this context, remote sensing-based methods have already been proven as an effective alternative for mapping rice area [6–9] and land use/land cover [10–12] as well.

During the last two decades, remote sensing-based methods have been widely used in mapping the rice cultivated areas. The most commonly used methods are developed on the basis of exploiting the optical remote sensing-based spectral indices. The implementation of these indices has been observed in many countries and some of the example cases are summarized in

Table 1. Examples of vegetation indices appearing in the literature used for mapping rice areas.

Ref.	Approaches
[13]	Utilized SPOT VGT-derived NDVI images over Mekong delta, Vietnam. They employed unsupervised classification and divergence statistics to determine signature separabilities; and found overall accuracy of 94%.
[4,14]	Employed MODIS-derived LSWI, NDVI, and EVI images over Southern China, and South and Southeast Asia. They introduced rice mapping algorithm based on the relationships between LSWI, NDVI and EVI during the plantation time; and obtained reasonable agreement with Landsat ETM+ derived rice maps [14], and national agricultural statistical data [4] (r2 ranging from 0.42–0.87).
[15]	Used MODIS-derived NDVI-values over six South Asian countries (i.e., India, Pakistan, Nepal, Bangladesh, and Bhutan). They applied several methods, such as spectral matching, decision trees, and perfect temporal profile associated with rice growth stages; and obtained strong relations with agricultural census data (i.e., r2 ≥ 0.97).
[16]	Utilized MODIS-derived LSWI and EVI images over Hunan, China for mapping rice cropping patterns that include single-season rice, early-season rice, and late-season rice systems. They applied a variable EVI/LSWI threshold function to recognize rice fields at the flooding stage. Outcomes were validated with finer spatial resolution SPOT-5 HRG (i.e., 2.5 m) land use information, and agriculture data on the country-level; and showed fair results (i.e., r2 of 0.52 for early season rice, 0.59 late season rice, and 0.34 for single season rice).
[17]	Used MODIS-derived NDVI, RVI and SAVI over Bali, Indonesia. They applied temporal variance analysis of the vegetation indices to clusterize between rice areas from other land uses; and obtained high correlation with agricultural data (i.e., r2 of 0.97 for regional-level and 0.92 for district level).
[18]	Used MODIS-derived NDVI, EVI, LSWI, and NDBI over Indonesia. They developed an algorithm based on temporal profile of vegetation strength and water content; and obtained a strong agreement with national-level agricultural data (i.e., r2 of 0.97).
[19]	Employed Landsat ETM+ over Bali, Indonesia. They evaluated the relationship between rice age, rice spectral and vegetation indices (that include: NDVI, RVI, IPVI, DVI, TVI, SAVI, and RGVI) and demonstrated significant relations with agricultural data (i.e., r2 of 0.97).
[20]	Used AVHRR-derived SAVI-values and Landsat TM-derived rice area over Hubei, China. They applied maximum likelihood classifier and found strong relations with agricultural data (i.e., an overall accuracy in the range 84.5–91.6%).

Note: AVHRR—Advanced Very High Resolution Radiometer; DVI—difference vegetation index; EVI—enhanced vegetation index; IPVI-infrared percentage vegetation index; LSWI—land surface water index; MODIS—Moderate Resolution Imaging Spectroradiometer; NDBI—normalized build-up index; NDVI—normalized difference vegetation index; NOAA—National Oceanic and Atmospheric Administration; RGVI—rice growth vegetation index; RVI—ratio vegetation index; SAVI—soil adjusted vegetation index; SPOT-5 HRG—Satellite Pour l’Observation de la Terre 5 High Resolution Geometric; SPOT VGT—Satellite Pour l’Observation de la Terre (SPOT) vegetation; TVI- transformed vegetation index.

.

In most of the abovementioned studies (see Table 1 [4,13–20]), the spatial resolution of the input images (i.e., MODIS, AVHRR, and SPOT-VGT) were having relatively low (i.e., in the range 500 m to 1.1 km) [4,13–18]. This particular issue would be critical as the dimension of some rice fields might be smaller than those of the employed spatial resolution. In this context, the use of Landsat imagery would compensate the spatial resolution [19]; however, it would be very difficult to obtain cloud-free chronology of such images over the entire growing season. In this paper, our overall objective was to develop a methodology for mapping boro rice (i.e., cultivated during the months January to May) using MODIS-derived NDVI values at 250 m spatial resolution and its implementation over Bangladesh during the 2007–2012 period. The specific objectives included the: (i) Determination of rice signatures in temporal dimension; (ii) formulation of a mathematical model for extracting rice cultivated areas; (iii) model calibration and validation; and (iv) delineation of the spatial distribution of boro rice.

2. Materials

2.1. General Description of the Study Area

Our study area is the entire country of Bangladesh, which is situated in the Indo-Gangetic plain of Southeast Asia. It is located on the northern edge of the Bay of Bengal, bordered by India on three sides (i.e., north, west, and a portion of the east), and also partly by Myanmar on the east side (Figure 1a). Geographically, it lies between 20°34′ to 26°38′N Latitude and 88°01′ to 92°41′E Longitude (Figure 1b).

Topographically, it is mostly flat 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. Climatically, the study area experiences sub-tropical climatic conditions with three prominent seasons: the dry winter (December to February), the pre-monsoon hot summer (March to May), and the rainy monsoon season (June to October) [21]. Annually, January is the coldest month (i.e., temperature is ∼17 °C) and April is the hottest one (i.e., temperature ranging from 30–40 °C). Bangladesh receives a minimum of 2000 mm of rainfall a year for most of the area except the northwestern region of Rajshahi (i.e., ∼1600 mm), and coastal area of Chittagong in the southeast and northeastern part of Sylhet (i.e., ∼4000 mm) [22]. Note that ∼80% of the rainfall occurs during the monsoon season (i.e., June to October) and very little during the dry season (i.e., December to February). The country experiences natural calamities, such as floods due to heavy monsoon rainfall, droughts during dry season and tropical cyclones originating over the Bay of Bengal during the period April to May and October to November. All of these disasters are prime causes of significant losses of agricultural productivity.

The agriculture sector in Bangladesh (i.e., accounts for about 62.8% of total area) contributes about 35% of the Gross Domestic Product (GDP) and employs about 70% of the labor force [23]. About 70% of agricultural land is dedicated to rice production, which constitutes approximately 92% of total food grains produced annually and accounts for nearly 18% of the GDP and provides about 60–70% of total calorie intake [24,25]. Bangladesh has three types of rice crops/seasons within a year, namely: (i) Aman with the growing season between June/July and November/December that depends on the rainfall until mid-September and then on the irrigation; (ii) aus between March/April and July/August that grows mainly under the influence of rainfall; and (iii) boro between January/February and May/June that primarily grows under irrigation scheme. Note that boro is the highest contributors of total rice supply in Bangladesh (i.e., about 55% of the total rice production during July 2011 to June 2012 period [26]).

2.2. Data Requirements and Its Pre-Processing

In this study, we primarily used MODIS-derived 16-day composite of NDVI (i.e., MOD13Q1 version 005) at 250 m spatial resolution, which were rectified through geometric correction, atmospheric calibration, and validated and quality assured through the earth observing system program [27,28]. These images were acquired during the boro growing season (i.e., January to May) for the years 2007–2012; and required some level of pre-processing prior to the development of the rice mapping system. Those included: (i) Re-projection of the NDVI images into Transverse Mercator from the default Sinusoidal Gridded Projection system; (ii) exclusion of cloud-contaminated pixels using Quality Assessment information; (iii) mosaic of the two adjacent NDVI images in order to have continuous images for the entire study area; and (iv) generation of the NDVI time-series for each of the years.

In addition, we used two sources of information to aid in the identification of boro cropping areas, such as (i) very high spatial resolution imagery with <5 m spatial resolution available from Google Earth during 2010, which was found to be effective for the visual interpretation of land cover and identification of the rice crop fields in other studies [15]; and (ii) agricultural land use map available from Bangladesh Soil Resource Development Institute. We also acquired ground-based boro rice cultivation area information from Bangladesh Bureau of Statistics (BBS) during the period 2007–2012 at 23 district levels (see Figure 1b for their spatial extents). These data were collected through field visits and interviewing the farmers, and then aggregated at 23 district levels; and used to calibrate and validate the MODIS-based boro rice mapping system.

3. Methods

Figure 2 shows a schematic diagram of the proposed method for developing a boro rice mapping system. It consisted of three major components, (i) ISODATA clustering and determining the boro rice signatures in temporal dimension; (ii) formulating a mathematical model for extracting boro rice areas; and (iii) model calibration and validation. These are briefly described in the following subsections.

3.1. ISODATA Clustering and Determining the Boro Rice Signatures

We determined the NDVI-based temporal signatures related to boro rice growing season in two steps. In the first step, we applied the ISODATA clustering technique over the NDVI time-series for each of the seasons for the period 2007–2009 in order to determine boro rice-specific signatures. For this technique, we considered having 40 classes by setting the convergence threshold of 0.995 (i.e., similar to other studies [9,15]) with the infinite number of iterations (i.e., the iteration would stop upon reaching to the required convergence threshold). Note that the rationale behind the adoption of this particular clustering was its ability to: (i) Assign every pixel to a class based on spectral and temporal similarities [29]; (ii) define many classes based on inherent characteristics of the time-series; (iii) classify a large area covered by various vegetation types (i.e., heterogeneous landscape) [15].

In the next step, the temporal signature (i.e., mean NDVI-values over the entire cropping season) for each of the classes was individually evaluated against various ancillary data. These included: (i) Land use map to make sure the class of interest fell under agricultural land category; (ii) very high spatial resolution imagery to determine visually identifiable unique features related to rice plantation, such as small dikes known as bunds, irrigation structures (that include canals, irrigation channels, ponds, etc.) [15]; and (iii) boro rice phenological information [30] to find agreement with the NDVI temporal dynamics for the class of interest (i.e., NDVI-values would be low at the start and end of the growing season; and high at the peak of the season). The signatures satisfying the abovementioned conditions would eventually be defined as the generic set of “boro rice signatures”.

3.2. Formulating Mathematical Model for Extracting the Boro Rice Areas

We formulated a mathematical model based on exploiting the generic set of signatures generated in Section 3.1 (see Figure 2b for more details). For each of the classes we computed signature-specific average A_i and standard deviation SD_i for the entire NDVI time-series. Then, we considered all of the signatures and deduced the following:

$$A=\frac{\sum _{i=1}^n{A}_i}{n}$$

(1)

$${\mathit{SD}}_A=\sqrt{\frac{1}{n-1}\sum _{i=1}^n{\left(A_i-A\right)}^2}$$

(2)

$$A_{\mathit{SD}}=\frac{\sum _{i=1}^n{\mathit{SD}}_i}{n}$$

(3)

$${\mathit{SD}}_{\mathit{SD}}=\sqrt{\frac{1}{n-1}\sum _{i=1}{\left({\mathit{SD}}_i-A_{\mathit{SD}}\right)}^2}$$

(4)

where, A and SD_A are the average and standard deviation respectively for all of the class-specific A_i-values;

A_SD and SD_SD are the average and standard deviation respectively for all of the class-specific SD_i-values; and n is the total number of classes.

In the course of implementing the abovementioned mathematical model, we computed both average and standard deviation images for each of the years using a 16-day composite of NDVI images. In these images, we denoted each of the pixel/cells as P_avg and P_std for the average and standard deviation images, respectively. Finally, we determined a pixel to be labeled as boro rice cultivation, if it satisfied the following set of formula. Note that these were similar to that of Nuarsa et al. [17].

$$(A-M*{\mathit{SD}}_A)\le P_{\mathit{avg}}\le (A+M*{\mathit{SD}}_A)$$

(5)

$$(A_{\mathit{SD}}-M*{\mathit{SD}}_{\mathit{SD}})\le P_{\mathit{std}}\le (A_{\mathit{SD}}+M*{\mathit{SD}}_{\mathit{SD}})$$

(6)

where, M is a multiplication factor, and its optimal magnitude can be determined through the calibration process described in the following sub-section.

3.3. Model Calibration and Validation

In calibration phase, we determined the optimal value of M to be used in Equations (5) and (6). In this case, we used the NDVI time-series derived year-specific average and standard deviation images during the period 2007–2009 (i.e., 50% of the data). For each of the years, we changed the values of M in the range of 1–1.5, which would magnify both standard deviations (i.e., SD_A and SD_SD) in Equations (5) and (6). The optimal value of M was determined when we observed the best agreements between the model (i.e., MODIS-based) and ground-based estimates of boro rice area at the country-level. Also, we performed linear regression analysis at 23 district levels to evaluate the effectiveness of the optimal M value.

In validation phase, we applied the optimal value of M (as determined in the last paragraph) over an independent dataset; i.e., NDVI time-series derived year-specific average and standard deviation images during the period 2010–2012 (i.e., 50% of the data). This process led to producing the spatial extent of the boro rice area for each of the individual years during the period 2010–2012. Then, we calculated the boro rice area over 23 districts and compared with the ground-based information. The degree of agreements between the model and ground-based estimates of boro rice area were determined in two ways, by means of calculating: (i) percentage error at country-level; and (ii) coefficient of determination (i.e., r²) upon performing linear regression at 23 district levels.

4. Results and Discussion

4.1. Determining the Boro Rice Signatures

Figure 3 shows the extracted signatures during the period 2007–2009. For each of the years, we found that 12 classes in 2007 (see Figure 3a), 13 in both 2008 and 2009 (see Figure 3b,c respectively), illustrated the boro rice area. In general, the signatures from the various classes that represented boro rice were similar to each other in all of the years. They also demonstrated distinct characteristics during three boro crop development stages. These included: (i) The initial/transplanting stage spanning between day of year (DOY) 01–32 (i.e., 1 January to 1 February) showed the lowest NDVI-values as the amount of biomass was at a very low level; (ii) the peak stage from mid-March to mid-April (that correspond to DOY from 81–112 to a large extent) depicted the highest NVDI-values due to the fact that the crop reached the maximum level of greenness; and (iii) the harvesting stage took place from DOY 129–161 (i.e., 9 May to 9 June) and again demonstrated very low NDVI-values. Between the initial and peak stages, the NDVI-values showed an increasing trend, due to the gradual development of the crop towards its maximum greenness/maturity. On the contrary, the NDVI-values showed a decreasing trend between the peak and harvesting periods as crops initiated the process of ripening (i.e., changing its color from green to yellowish/goldish with little to no chlorophyll content). Our observed developmental stages were similar to those reported in the literature (e.g., [30]). Also, we found that remote sensing-derived signatures over rice crops in other countries demonstrated similar developmental stages [14,16,31–34].

4.2. Model Calibration and Validation

During the calibration phase, we found that the optimal value of M was 1.25, which revealed strong agreements between the model and ground-based estimates of boro rice area at the country-level during the period 2007–2009 (see

Table 2. Determination of optimal value of multiplier illustrated in Equations (5) and (6) during 2007–2009 at the country-level boro rice cultivation area.

Year	Ground-Based Estimates (ha)	MODIS-Based Estimates (ha)
		Value of the Multiplier, M
		1	1.1	1.2	1.25	1.3	1.4	1.5
2007	4,257,873	2,816,793	3,360,706	3,936,150	4,237,718	4,535,018	5,153,543	5,772,318
2008	4,607,630	2,832,081	3,359,081	3,914,050	4,493,881	4,787,100	5,086,418	5,692,337
2009	4,716,247	3,449,918	4,061,175	4,313,650	4,695,150	5,348,693	6,012,843	6,670,575

). The application of this M-value during the 2007–2009 period also demonstrated strong linear relationships (i.e., r² in the range 0.84–0.91) compared to the ground-based estimates at the 23 district levels (see Figure 4).

During the validation phase, we implemented the optimal M-value (= 1.25) over an independent dataset during 2010–2012; and found that the level of agreements was strong (i.e., relative error in the range −0.83–1.42%) between the model and ground-based estimates at the country-level (see

Table 3. Application of the optimal value of the multiplier illustrated in Equations (5) and (6) during 2010–2012; and comparison between MODIS and ground-based estimates boro rice cultivation area at the country level.

Year	Ground-Based Estimates (ha)	MODIS-Based Estimates (ha)	Relative Error
2010	4,706,875	4,639,975	1.42
2011	4,770,337	4,757,018	0.28
2012	4,810,025	4,850,062	−0.83

). In addition, Figure 5 illustrates linear regression analysis between the model and ground-based estimates of boro rice at 23 district levels during the period 2010–2012. It revealed strong relations (i.e., r² in the range of 0.69–0.89) between the variables of interest. Note that similar strong relations between remote sensing and ground-based estimates were also reported in other studies over different countries, such as (i) Xiao et al. [4] who found r² to be in between 0.42 and 0.87 in South and Southeast Asia; (ii) Gumma et al. [15] demonstrated strong relationships (i.e., r² ≥ 0.97) over six South Asian countries; (iii) Fang et al. [20] observed strong relations (i.e., an overall accuracy of 84.5–91.6%) over China; and (iv) Bridhikitti and Overcamp [34] found r² to be in between 0.66 and 0.93 over Southeast Asia. In addition, we also found the root mean square errors were (i) 71,763 ha in 2010; (ii) 43,844 ha in 2011; and (iii) 48,925 ha in 2012. In addition, the calculation of mean relative errors was: (i) 0.40% in 2010; (ii) 0.44% in 2011; and (iii) −1.36% in 2012 over 19 districts except the districts of Bandarban, Khagrachhari, Patuakhali and Rangamati. In these four districts, we found a mixing of the boro rice signatures with other land use/land cover types, which might require further investigation using higher spatial resolution images.

In both calibration and validation phases, we observed that the correlation between the model and ground-based estimates of boro area was captured to a large extent (i.e., r² in the range of 0.75–0.91). We assumed that unexplained differences were associated with one or more of the following factors. For example: (i) The spatial resolution (i.e., 250 m) of the employed MODIS data might not be suitable for all of the crop fields as some of them were smaller and also heterogeneous in nature; (ii) the initial/transplanting stage was assumed to be the month of January, however it might be possible to have a week or two shifting in some places; and (iii) the ground-based estimates of the rice area were collected through field visits and interviewing the farmers; thus some error might be induced during this process.

4.3. Spatial Distribution of Boro Rice

Figure 6 shows an example spatial distribution of boro rice during 2012. Usually, the spatial distribution of this crop depends on several interrelated factors, such as climate (i.e., temperature, humidity), soil conditions (i.e., nutrient availability, soil texture, soil structure, salinity), topography, and management practices (i.e., water supply, fertilizer use) [35,36].

In general, boro rice was produced very well (i.e., at least 50% of the areas under boro cultivation) in all of the districts except for six, and is summarized as follows:

• Rajshahi district in the northwest produced less due to the existence of Barind Tract (i.e., consisting of terraced lands and is also the driest region in the country) [37,38];
• Khulna in the southwest observed less production, where the dominant land use/land cover was mangrove forest and shrimp cultivation [38,39]; and
• Khagrachhari, Rangamati, and Bandarban in the southeast; and Sylhet in the northeast also produced less due to the mountain ranges [38,39].

5. Conclusions

In this paper, we demonstrated a simple protocol (that included ISODATA clustering technique; formulation of the mathematical model; and model calibration and validation) for mapping the boro rice area and its implementation over Bangladesh using a MODIS-derived 16-day composite of NDVI. However, we would only require the mathematical model (i.e., system) in rice mapping using the NDVI time-series directly, which would produce the map within a couple of minutes. This particular aspect of our proposed system could be considered as a unique advancement over other methods. In the context of the outcomes, our analysis showed that the developed system provided a reasonable agreement with ground-based estimates at the country and district levels in the period from 2007–2012, both in spatial distribution and quantity. Regression analysis showed r² values ranged from 0.84–0.91 and from 0.75–0.89 during the period of 2007–2009 and 2010–2012, respectively, at 23 district levels. In addition, we believe that the proposed model and MODIS data used in this study would be useful in: (i) Mapping the rice areas at large scale; (ii) generating an updated and timely database for rice cultivated areas, which are crucial to research on rice production impacts on methane emissions and quality and amount of water resources; and (iii) mapping rice areas, where the ground truth data are not available. Despite the effectiveness of the employed method, we suggest that the proposed model should be calibrated and validated prior to implementation in other regions across the world. In addition, it would be worthwhile mentioning that our proposed method would be extendable to those areas with a small proportion of rice in a cultivated area, if the rice fields are detectable in 250 m spatial resolution images and the phenology of the co-existing crops are different from the rice.