An Object-Based Paddy Rice Classification Using Multi-Spectral Data and Crop Phenology in Assam , Northeast India

Rice is the staple food for half of the world’s population. Therefore, accurate information of rice area is vital for food security. This study investigates the effect of phenology for rice mapping using an object-based image analysis (OBIA) approach. Crop phenology is combined with high spatial resolution multispectral data to accurately classify the rice. Phenology was used to capture the seasonal dynamics of the crops, while multispectral data provided the spatial variation patterns. Phenology was extracted from MODIS NDVI time series, and the distribution of rice was mapped from China’s Environmental Satellite (HJ-1A/B) data. Classification results were evaluated by a confusion matrix using 100 sample points. The overall accuracy of the resulting map of rice area generated by both spectral and phenology is 93%. The results indicate that the use of phenology improved the overall classification accuracy from 2%–4%. The comparison between the estimated rice areas and the State’s statistics shows underestimated values with a percentage difference of  ́34.53%. The results highlight the potential of the combined use of crop phenology and multispectral satellite data for accurate rice classification in a large area.


Introduction
Studies on the extent of rice areas provide useful information for food security, water resources management and environmental sustainability.Rice is the major food for nearly half of the world's seven billion people, mostly in developing countries in Asia, Africa and Latin America [1].Rice agriculture is mostly irrigated and consumes 24%-30% of world's developed fresh water resources [2].Therefore, determining the total rice area is an important input for the effective management of global water resources.On the other hand, the world's rice production is not increasing significantly, and present annual rice demand exceeds annual production; thus, food security remains uncertain [3].Rapid urbanization, industrialization, changing patterns of precipitation and rising global temperature affect the land and water resources for rice production [4].Hence, it is urgent to monitor the rice area to meet the growing food demand, efficient water resource management and environmental sustainability.
Remote sensing (RS) has demonstrated its potential for rice area mapping employing either optical or synthetic aperture radar (SAR) images [5][6][7][8][9][10].Time series multi-spectral and multi-temporal data from MODIS, NOAA-AVHRR, LANDSAT, IRS, SPOT and Chinese HJ-1A/B were used by many researchers at the national and sub-national scale for rice mapping [11][12][13][14][15]. Though optical data cover a wide range of spatial and temporal resolutions, it is weather dependent, and cloud cover hampers its operation.To overcome this drawback, weather-independent SAR has been used [9,10,16].
Unfortunately, the low temporal resolution and high price of SAR data limit its use at a large scale.A number of techniques have been applied to discriminate rice area, mainly the thresholding method [17][18][19], the supervised classification method [20,21] and phenology-based mapping [11,22].However, these methods have limitations.Thresholding requires appropriate thresholding values for accurate classification.Supervised classification requires training data for each year, which is a challenge for large area mapping.The phenology-based mapping requires a long and continuous time series.Moreover, all of these methods underperform in the areas where the crop fields are small.
Compared to traditional pixel-based classification, object-based image analysis (OBIA) considers a group of pixels instead of a single pixel for classification [23].The OBIA produces more meaningful and reliable classification by contributing additional information, including spectral, textural and geometric features [24].The OBIA includes two steps: image segmentation and classification.In segmentation, the entire image is partitioned into regions or objects that are more internally uniform and homogeneous than neighboring objects [25].The segmented image generates extra spectral, textural and geometric information of objects.In the classification process, each object is assigned to a specific class according to its spectral, textural, geometric and customized properties.Previous studies have demonstrated the successful use of OBIA in many applications [24,[26][27][28].Small fragmented rice fields are common in India [29]; therefore, the OBIA approach is particularly useful by focusing on objects or parcels instead of concentrating only on the properties of single pixels.
Rice mapping for a large area is generally performed at low spatial resolution due to the insufficient availability of high resolution temporal data [17,19].To obtain the temporal data at high spatial resolution, image fusion or blending algorithms are commonly applied combing high and low spatial resolution images.The fusion algorithm can predict accurate surface reflectance [30] and has shown the potential of dense time series generation for phenology studies [31].
Phenology refers to the timing of growing events of plants, such as bud-burst, leafing, peak growth stage, flowering and abscission [32].In recent studies, phenology has been used for crop classification [11,22,33,34].However, the capability of phenology has not been investigated yet for the OBIA-based rice classification.This study focuses on the use of phenology for rice classification under the OBIA framework.Its main objective is to investigate the applicability of phenology for OBIA-based rice mapping.

Study Area
The study area includes five districts (Bongaigaon, Barpeta, Goalpara, Nalbari and Kamrup) of Assam in northeast India, with an area about 14,000 square km.The center of the area is located at 26 ˝23 1 N and 91 ˝09 1 E (Figure 1).Assam is mainly the flood plain of two rivers: Brahmaputra and Barak.The region experiences heavy rainfall with an average of 3000 mm annually.Rice, the principal crop of the region, is cultivated during very wet summers, as well as in very dry winters.The area under rice cultivation is 25 million hectares, which is 71% of the total cultivated area [35].
The crop is cultivated in many environmental conditions in the region based on the hydrological characteristics.The most commonly-practiced cultivation types are: (1) rainfed lowland rice, where the flooding of the fields is non-continuous and water depth varies around 1 cm-50 cm; the water availability of the fields primarily depends on the rainfall; (3) irrigated rice, where the water depth of the fields varies between 5 cm and 10 cm; the fields are kept constantly flooded throughout the growing seasons; two crops are harvested per year from these fields, and often, the production is high; (3) flood-prone rice, which is cultivated in deep-water areas with a water depth of 0.5 m-3 m; this rice is generally grown near the rivers and submerged frequently; the production is relatively low due to the effect of floods; and (4) upland rice, which is mostly cultivated on the hill slopes of the region; it is practiced in the rainfed conditions without the flooding and accumulation of water in the fields; the water requirement is low for this type of rice.
The region has two distinct growing seasons of rice: rabi (February/March-June/July) and kharif (June/July-November/December) (Figure 2).Most farmers prefer to plant kharif rice due to the abundant rain and favorable temperature.Rabi rice is mainly distributed near the river and the lakes to take advantage of the captured or retained water of the monsoon rain or the irrigation infrastructure.The distribution of double-cropped rice is scattered on the north bank of Brahmaputra River, where irrigation is available.The selected study area has all of the rice types, and all were considered for the mapping.to take advantage of the captured or retained water of the monsoon rain or the irrigation infrastructure.The distribution of double-cropped rice is scattered on the north bank of Brahmaputra River, where irrigation is available.The selected study area has all of the rice types, and all were considered for the mapping.

MODIS Data
Moderate Resolution Imaging Spectroradiometer (MODIS) 16-day composited Normalized Difference Vegetation Index (NDVI) products (MOD13Q1, Collection 5) were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) [36].MODIS NDVI is a 250-m spatial resolution 16-day composite dataset.The MODIS NDVI is suitable for crop monitoring with 23 images per year.All of the available MODIS NDVI data were acquired from April 2014-March 2015.All of the datasets were re-projected to WGS 84/UTM Zone 36 N and then resampled to 30-m spatial resolution to make them compatible with the HJ-1A/B datasets.

HJ-1A/B Data
The Chinese HJ-1 A/B satellite's CCD (charge coupled device) sensors' data acquired during the growth period of rice for the year of 2014-2015 were used in this study (Table 1).The HJ-1 A/B CCD images have good quality and can be used to extract crop areas [5,37].This dataset is available for download from CRESDA (China Centre for Resources Satellite Data and Application) [38].The sensors characteristics of this satellite are given in Table 2.The radiometric calibration was performed using the coefficients provided by the CRESDA [38].The atmospheric correction of the images was performed using the FLAASH atmospheric correction model of ENVI 5 [39].to take advantage of the captured or retained water of the monsoon rain or the irrigation infrastructure.The distribution of double-cropped rice is scattered on the north bank of Brahmaputra River, where irrigation is available.The selected study area has all of the rice types, and all were considered for the mapping.

MODIS Data
Moderate Resolution Imaging Spectroradiometer (MODIS) 16-day composited Normalized Difference Vegetation Index (NDVI) products (MOD13Q1, Collection 5) were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) [36].MODIS NDVI is a 250-m spatial resolution 16-day composite dataset.The MODIS NDVI is suitable for crop monitoring with 23 images per year.All of the available MODIS NDVI data were acquired from April 2014-March 2015.All of the datasets were re-projected to WGS 84/UTM Zone 36 N and then resampled to 30-m spatial resolution to make them compatible with the HJ-1A/B datasets.

HJ-1A/B Data
The Chinese HJ-1 A/B satellite's CCD (charge coupled device) sensors' data acquired during the growth period of rice for the year of 2014-2015 were used in this study (Table 1).The HJ-1 A/B CCD images have good quality and can be used to extract crop areas [5,37].This dataset is available for download from CRESDA (China Centre for Resources Satellite Data and Application) [38].The sensors characteristics of this satellite are given in Table 2.The radiometric calibration was performed using the coefficients provided by the CRESDA [38].The atmospheric correction of the images was performed using the FLAASH atmospheric correction model of ENVI 5 [39].

MODIS Data
Moderate Resolution Imaging Spectroradiometer (MODIS) 16-day composited Normalized Difference Vegetation Index (NDVI) products (MOD13Q1, Collection 5) were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) [36].MODIS NDVI is a 250-m spatial resolution 16-day composite dataset.The MODIS NDVI is suitable for crop monitoring with 23 images per year.All of the available MODIS NDVI data were acquired from April 2014-March 2015.All of the datasets were re-projected to WGS 84/UTM Zone 36 N and then resampled to 30-m spatial resolution to make them compatible with the HJ-1A/B datasets.

HJ-1A/B Data
The Chinese HJ-1 A/B satellite's CCD (charge coupled device) sensors' data acquired during the growth period of rice for the year of 2014-2015 were used in this study (Table 1).The HJ-1 A/B CCD images have good quality and can be used to extract crop areas [5,37].This dataset is available for download from CRESDA (China Centre for Resources Satellite Data and Application) [38].The sensors characteristics of this satellite are given in Table 2.The radiometric calibration was performed using the coefficients provided by the CRESDA [38].The atmospheric correction of the images was performed using the FLAASH atmospheric correction model of ENVI 5 [39].

Field Survey Data
In May 2015, a field survey was conducted by using the GVG (GPS-video-GIS) instrument developed by Wu et al. [40].GVG is capable of producing the crop area estimation and crop type proportion of the surveyed area with very high accuracy [41].A distance of 400 km was covered across the study area during the survey, and around 1800 ground truth points were collected (Figure 1).The ground truth points include the rice fields and other major land cover classes.

Agriculture Statistics Data
Agricultural statistics at the district level were obtained from the Directorate of Economics and Statistics, Government of Assam [42].The dataset contains district-level estimation of cropped area, yield and production for the different crops.The data of the rice areas were used to check the agreement between the remote sensing-based results and the government rice area statistics.

Ground Reference Data
The following auxiliary information has been used to generate the ground reference map: (1) LULC (land use/land cover) map of 2013-2014 (scale: 1:250 K); (2) LULC: SIS-DP (space-based information support for decentralized planning) (scale: 1:10 K); (3) field survey data; and (4) Google Earth images.Datasets (1) and ( 2) are available at Bhuvan [43].The datasets were integrated into the systems as OGC Web Services using the open source geographic information system QGIS [44].The Google Earth images helped to identify crop types.The final ground reference map was created after the analysis of above-mentioned data and manually digitizing all of the homogeneous rice fields for the years of 2014-2015.The dataset was utilized in training sample selection and accuracy assessment.

Methods
The flowchart of the method is shown in Figure 3.The major steps include data pre-processing, fusion of the MODIS and HJ NDVI, phenology extraction from the NDVI time series, segmentation of the images, classification, as well as validation.The main feature of the method is that the phenological variables were extracted from the fused time series high resolution NDVI data and combined with the HJ CCD spectral data for OBIA-based rice classification.The approach utilized not only the seasonal features through phenological variables, but also the spectral features of HJ CCD data.

Construction of Smooth Time Series NDVI
Although MODIS NDVI is calculated from atmospherically-corrected surface reflectance and has reduced noise, cloud and aerosol effects, the time series exhibits noise and produces inaccurate phenology.The Savitzky-Golay (S-G) filter [45] was applied to smooth the data.This is achieved by applying a locally-adaptive moving window, which uses a polynomial least square regression to fit the data.In Figure 4, smoothed NDVI is shown for an agriculture landscape.

Construction of Smooth Time Series NDVI
Although MODIS NDVI is calculated from atmospherically-corrected surface reflectance and has reduced noise, cloud and aerosol effects, the time series exhibits noise and produces inaccurate phenology.The Savitzky-Golay (S-G) filter [45] was applied to smooth the data.This is achieved by applying a locally-adaptive moving window, which uses a polynomial least square regression to fit the data.In Figure 4, smoothed NDVI is shown for an agriculture landscape.

Construction of Smooth Time Series NDVI
Although MODIS NDVI is calculated from atmospherically-corrected surface reflectance and has reduced noise, cloud and aerosol effects, the time series exhibits noise and produces inaccurate phenology.The Savitzky-Golay (S-G) filter [45] was applied to smooth the data.This is achieved by applying a locally-adaptive moving window, which uses a polynomial least square regression to fit the data.In Figure 4, smoothed NDVI is shown for an agriculture landscape.

Blending of Time Series MODIS NDVI and HJ NDVI
An enhanced spatial and temporal adaptive reflectance fusion (ESTARFM) [30] was selected to fuse MODIS and HJ NDVI.The algorithm uses two or more pairs of fine and coarse resolution images acquired on the same date and a coarse resolution image acquired on the prediction date.It predicts the pixel values based on the temporal weight determined by the spectral similarity between fine and coarse pixels [30].To obtain an improved prediction, three HJ NDVI images acquired at key rice growth stages covering the whole growing seasons were utilized.The HJ NDVI images were acquired 22 October 2014, 5 December 2014 and 9 March 2015, which corresponds to the flowering, matured ripened and heading stages of rice, respectively.The HJ NDVI pairs nearest to the MODIS NDVI date were selected for the fusion.The index-then-blend (IB) approach was adopted for the fusion, where vegetation indices were calculated first and then using these indices for the fusion [46].The fusion produced time series MODIS NDVI at enhanced spatial resolution of 30 m. Figure 5 shows an example of a MODIS input image and a corresponding synthetic MODIS image.The synthetic image will be similar, like the nearest time 30-m HJ CCD NDVI image with very little expected change in 5 days.The ESTARFM is computationally intensive [30]; for this study, it took around 103 computing hours with an Intel Core i7, 16 GB RAM, Windows 7 Computer.Therefore, to save computation time, it is suggested to use parallel processing or GPU (graphics processing unit)-enabled processing.

Blending of Time Series MODIS NDVI and HJ NDVI
An enhanced spatial and temporal adaptive reflectance fusion (ESTARFM) [30] was selected to fuse MODIS and HJ NDVI.The algorithm uses two or more pairs of fine and coarse resolution images acquired on the same date and a coarse resolution image acquired on the prediction date.It predicts the pixel values based on the temporal weight determined by the spectral similarity between fine and coarse pixels [30].To obtain an improved prediction, three HJ NDVI images acquired at key rice growth stages covering the whole growing seasons were utilized.The HJ NDVI images were acquired 22 October 2014, 5 December 2014 and 9 March 2015, which corresponds to the flowering, matured ripened and heading stages of rice, respectively.The HJ NDVI pairs nearest to the MODIS NDVI date were selected for the fusion.The index-then-blend (IB) approach was adopted for the fusion, where vegetation indices were calculated first and then using these indices for the fusion [46].The fusion produced time series MODIS NDVI at enhanced spatial resolution of 30 m. Figure 5 shows an example of a MODIS input image and a corresponding synthetic MODIS image.The synthetic image will be similar, like the nearest time 30-m HJ CCD NDVI image with very little expected change in 5 days.The ESTARFM is computationally intensive [30]; for this study, it took around 103 computing hours with an Intel Core i7, 16 GB RAM, Windows 7 Computer.Therefore, to save computation time, it is suggested to use parallel processing or GPU (graphics processing unit)-enabled processing.

Image Segmentation and Object Feature Selection
The main purpose of this study is to identify area cropped with rice.Rice fields can be extracted as an object or segment by applying image segmentation methods.Image segmentation generates the segments as regions or groups of pixels with certain homogeneity criteria and includes additional spectral, spatial and textural information as compared to a single pixel [23].Our study utilized these additional features for rice crop classification.The detailed description of an object's features can be read in Definiens [47].However ,spatial and textural features have not been considered, as they do

Image Segmentation and Object Feature Selection
The main purpose of this study is to identify area cropped with rice.Rice fields can be extracted as an object or segment by applying image segmentation methods.Image segmentation generates the segments as regions or groups of pixels with certain homogeneity criteria and includes additional spectral, spatial and textural information as compared to a single pixel [23].Our study utilized these additional features for rice crop classification.The detailed description of an object's features can be read in Definiens [47].However ,spatial and textural features have not been considered, as they do not contribute much for agricultural crop identification [48]; textural features result in longer computation time and provide only minor improvement in classification results [24,28].The multi-resolution segmentation algorithm [49] implemented in eCognition software has been used to generate the rice field objects.Cloud-free HJ-1A/B images with spectral Bands 2, 3 and 4, i.e., green, red and near infra-red (NIR), have been selected for the segmentation.Three segmentation parameters, shape, compactness and scale, were selected as 0.1, 0.7 and 25, respectively, after a series of tests.The value of the shape parameter determines the contribution of spectral weight.The compactness criteria decide the closeness of image objects.The scale parameter determines the size of the image objects.To investigate the optimum segmentation results, a pixel-based empirical discrepancy method [50] has been applied.For this purpose, a reference segmented image has been generated randomly, by ground truth objects covering the rice fields and compared to the algorithm-generated segmented image.The obtained average accuracy for pixel-based discrepancy was 90.8%.Moreover, the visual assessment of segmented results indicated comparable matching with the actual field sizes.The result of image segmentation is shown in Figure 6.At the last step in this process, for each segmented object, the mean value of all phenological variables and NDVI were calculated and assigned to each object.
Remote Sens. 2016, 8, 479 7 of 20 not contribute much for agricultural crop identification [48]; textural features result in longer computation time and provide only minor improvement in classification results [24,28].The multi-resolution segmentation algorithm [49] implemented in eCognition software has been used to generate the rice field objects.Cloud-free HJ-1A/B images with spectral Bands 2, 3 and 4, i.e., green, red and near infra-red (NIR), have been selected for the segmentation.Three segmentation parameters, shape, compactness and scale, were selected as 0.1, 0.7 and 25, respectively, after a series of tests.The value of the shape parameter determines the contribution of spectral weight.The compactness criteria decide the closeness of image objects.The scale parameter determines the size of the image objects.To investigate the optimum segmentation results, a pixel-based empirical discrepancy method [50] has been applied.For this purpose, a reference segmented image has been generated randomly, by ground truth objects covering the rice fields and compared to the algorithm-generated segmented image.The obtained average accuracy for pixel-based discrepancy was 90.8%.Moreover, the visual assessment of segmented results indicated comparable matching with the actual field sizes.The result of image segmentation is shown in Figure 6.At the last step in this process, for each segmented object, the mean value of all phenological variables and NDVI were calculated and assigned to each object.

Derivation of Crop Phenology
The phenological variables were extracted from the fused 30-m NDVI time series using the TIMESAT toolbox [51], a widely-used software package for extracting seasonality from time series satellite data.The NDVI time series data were fitted to an asymmetric Gaussian (AG) function to extract the phenological variables [52].The parameter values of the AG model were adjusted interactively using the TIMESAT GUI to achieve the close fitting results.The seasonality parameter

Derivation of Crop Phenology
The phenological variables were extracted from the fused 30-m NDVI time series using the TIMESAT toolbox [51], a widely-used software package for extracting seasonality from time series satellite data.The NDVI time series data were fitted to an asymmetric Gaussian (AG) function to extract the phenological variables [52].The parameter values of the AG model were adjusted interactively using the TIMESAT GUI to achieve the close fitting results.The seasonality parameter value was selected as 0.1 in order to fit the double growing season of agricultural landscapes.To remove the spikes and outliers, the medium filter option with a parameter value of 2 was selected.To remove negatively-biased noise from NDVI, an upper envelope with a parameter value of 2 was selected, while adaptation strength was kept at the value of 1 to maintain the close fit result.
The phenology variables were extracted from NDVI time series by the commonly-used thresholding method, assuming that a specific phenology phenomenon occurs when NDVI values exceed a given threshold [52,53].The TIMESAT program iterates for each pixel of the time series, extracting the phenology variables according to the changing patterns of NDVI over time.In TIMESAT, the start and the end of the season is defined from the model function as the point in time for which the value has increased by a certain threshold [52].In this study, the threshold value is set to 10% of the distance between the minimum level and the maximum.
The large and small integrals (Table 3) are the annual integrated values of NDVI.These values are often used as a measure of net primary production [53,54].To achieve a good estimation of the vegetation production, only the growing season is considered for the integration of NDVI [52].
Among the other phenological variables extracted are the seasonal amplitude, the largest seasonal value, the seasonal base value and the length of the growing season.The meaning and the proxy information of extracted phenological variables is given in Table 3.In Table 3, only the first three lines describe actual phenological variables; the remaining ones are biomass related.In this study, all of the variables are referred to as a phenological variables.All of the listed variables in Table 3 are used in this study.
Table 3. Phenological variables and their definitions according to [52,53].The extracted phenology variables (Section 4.4) were composited with the HJ CCD spectral data.The composited data contained not only the spectral features of HJ CCD data, but also the phenological variables extracted from time series MODIS NDVI.A step-wise OBIA was performed on this composited dataset.An object-based image classification includes the process of image segmentation followed by assigning a specific class to the objects based on spectral, textural, geometric, as well as customized features.The entire OBIA classification algorithm was implemented in eCognition software [47].Under the OBIA framework, a top-down hierarchical approach has been adopted to identify the rice areas, as shown in Figure 7.At first, the broad general classes were identified and further sub-divided into smaller specific classes.The classification has been conducted mainly by applying an NN (nearest neighbor) classifier and a threshold to mean values of spectral features, NDVI and phenological features of image objects.In NN classification, at first, the classifier is trained by certain image objects as samples; then, all of the objects based on their nearest sample neighbors are classified [47].The threshold decides whether an image object qualifies as a specific condition or not for a specific class description [47].The classification of the images was carried out separately for the kharif and rabi seasons rice according to the crop calendar information.At first, the water and the land classes were distinguished by using mean NIR band values.The land was further sub-divided into vegetation and non-vegetation using NDVI.The vegetation was then classified into forest and non-forest by the mean value of the red band.Again, the non-forest was classified into cropland and grassland by the mean NDVI value.The grassland tends to have high NDVI all year with less seasonal changes.On the other hand, croplands have variable NDVI with high values only in the peak growing seasons.Further, the NN supervised classification was applied to sub-divide the cropland class into rice and non-rice (grassland, built-up and other crops) classes.The mean value of the red, green and NIR bands, as well as phenological variables were utilized to perform the NN classification.Finally, the classes other than rice were re-assigned to one "non-rice" class in order to obtain a final rice and non-rice classification.To analyze the effect of the phenology on rice classification, two classification approaches were used: (1) using only the multi-spectral bands;

Accuracy Assessment
To validate the thematic accuracy of object-based image classification, image objects should be used as sampling units [55], as this approach helps the evaluation of the entire classification process, including the segmentation [48].A confusion matrix was generated comparing classified image objects and reference image objects to evaluate the classification accuracy.In order to obtain the confusion matrix, a totally independent reference sample, i.e., 100 image objects (50 for rice; 50 for non-rice classes), was collected randomly with the help of ground reference data.Subsequently, overall accuracy, user's and producer's accuracy, as well as the kappa coefficient were calculated to assess the results [55,56].Additionally, derived rice area was compared to the agricultural statistics data.The percentage difference was measured to check the agreement between the estimated results and the rice area statistics.

Phenological Characteristics
Figure 8 shows the phenological variables derived from 16-day NDVI time series data according to the different vegetation types present in the region.The area under agriculture shows distinct

Accuracy Assessment
To validate the thematic accuracy of object-based image classification, image objects should be used as sampling units [55], as this approach helps the evaluation of the entire classification process, including the segmentation [48].A confusion matrix was generated comparing classified image objects and reference image objects to evaluate the classification accuracy.In order to obtain the confusion matrix, a totally independent reference sample, i.e., 100 image objects (50 for rice; 50 for non-rice classes), was collected randomly with the help of ground reference data.Subsequently, overall accuracy, user's and producer's accuracy, as well as the kappa coefficient were calculated to assess the results [55,56].Additionally, derived rice area was compared to the agricultural statistics data.The percentage difference was measured to check the agreement between the estimated results and the rice area statistics.

Phenological Characteristics
Figure 8 shows the phenological variables derived from 16-day NDVI time series data according to the different vegetation types present in the region.The area under agriculture shows distinct phenological patterns.The concept of phenology is pointless in non-vegetated areas, which are also characterized by low NDVI variability over time.The hills and evergreen vegetation area show a long growing season (~176 days) with an earlier start and a later end.The variable used as a proxy for the green biomass, such as base value (BV) and the largest NDVI value (LV), shows lower values in the cultivated land, whereas higher values area recorded in the high biomass regions.The large integral (LI) value, which is related to net primary production [53], shows the highest values in the cropped land and forested area as compared to other land cover classes.The seasonal amplitude (SA) and small integral (SI), which is linked to plant biomass [57], varied according to the forest cover type, as evident from the southern part of the region.All of the phenological variables demonstrated typical patterns for the various vegetation classes and showed potential for crop identification.However, further research is required to fully utilize all of the variables for crop classification.

Identified Rice with and without the Phenological Variables
The classification was performed for the two approaches: (1) using only the spectral feature and; (2) using both the phenological and spectral features, as shown in Figure 9.In each approach, rice fields were identified effectively.Classification performed with phenology increased the accuracy by 2%-3%.The confusion among grasslands, other crops and rice was minimized by the use of phenology.Using only one image, acquired at the ripening stage of rice (December), provided better results than the image of the heading stage (October).This is because mature ripened rice has typical spectral characteristics.The best classification result with an overall accuracy of 93% was achieved when all three images (October, December and March) were used.The changes of rice canopy structure over time is captured by using more images and increases the spectral separability.Figure 10 show the differences between the spectral and phenology + spectral classification results.It has been observed that the difference lies in the areas where vegetation has similar spectral characteristics, like rice crops.

Identified Rice with and without the Phenological Variables
The classification was performed for the two approaches: (1) using only the spectral feature and; (2) using both the phenological and spectral features, as shown in Figure 9.In each approach, rice fields were identified effectively.Classification performed with phenology increased the accuracy by 2%-3%.The confusion among grasslands, other crops and rice was minimized by the use of phenology.Using only one image, acquired at the ripening stage of rice (December), provided better results than the image of the heading stage (October).This is because mature ripened rice has typical spectral characteristics.The best classification result with an overall accuracy of 93% was achieved when all three images (October, December and March) were used.The changes of rice canopy structure over time is captured by using more images and increases the spectral separability.Figure 10 show the differences between the spectral and phenology + spectral classification results.It has been observed that the difference lies in the areas where vegetation has similar spectral characteristics, like rice crops.

Spatial Distribution of Rice
At 30-m spatial resolution, maps of the rice growing areas were generated (Figure 11).Most of the area cultivates kharif rice.The rabi rice is mainly distributed near the river and the lakes to take advantage of the captured or retained water of the monsoon rain or to benefit from irrigation facilities.However, rabi rice and double-cropped rice are less practiced in the region due to insufficient irrigation infrastructure.The distribution of double-cropped rice is scattered over the northern region.The map also shows that rice cultivated fields are small and medium in size over the region.

Spatial Distribution of Rice
At 30-m spatial resolution, maps of the rice growing areas were generated (Figure 11).Most of the area cultivates kharif rice.The rabi rice is mainly distributed near the river and the lakes to take advantage of the captured or retained water of the monsoon rain or to benefit from irrigation facilities.However, rabi rice and double-cropped rice are less practiced in the region due to insufficient irrigation infrastructure.The distribution of double-cropped rice is scattered over the northern region.The map also shows that rice cultivated fields are small and medium in size over the region.

Accuracy Assessment
The accuracy assessment showed a fairly high accuracy for the rice and the non-rice map (Table 4).Overall classification accuracy was above 84% for all of the cases when using phenological variables.The assessment showed that the overall accuracy was increased by 2%-4% by adding phenology with the spectral data.The kappa coefficient value was improved by 9%.The user accuracy for rice and non-rice was increased by 5% and 4.5%, respectively.The producer accuracy

Accuracy Assessment
The accuracy assessment showed a fairly high accuracy for the rice and the non-rice map (Table 4).Overall classification accuracy was above 84% for all of the cases when using phenological variables.The assessment showed that the overall accuracy was increased by 2%-4% by adding phenology with the spectral data.The kappa coefficient value was improved by 9%.The user accuracy for rice and non-rice was increased by 5% and 4.5%, respectively.The producer accuracy for both the rice and non-rice classes showed an improvement of approximately 4.3%.The errors are mainly due to the relatively small field sizes (less than one pixel).The rice fields are fragmented and usually mixed with other land cover classes, which leads to a mixed pixel problem and further decreases the classification accuracy.

Comparison with Agricultural Statistics Data
The estimated rice areas were compared to the State's rice area statistics at the district level (Table 5).The percentage difference was calculated for the estimated and the statistics rice areas to check the agreement between them.The estimated rice areas were comparable to the statistics for the Bongaigaon and Nalbari districts with a percentage difference of 9.31% and ´5.19%, respectively.However, derived rice areas were underestimated for the Barpeta, Goalpara and Kamrup with a percentage difference of ´53.68%, ´47.13% and ´51.29%, respectively.The discrepancy between the estimated and the statistics rice areas possibly is due to: (1) the limitations of the 30-m spatial resolution of HJ CCD in identifying small patches of paddy rice; and (2) where the compared datasets were from the different years and inconsistent statistics-generation methods for different districts.

Discussion
The integration of crop phenology has advantages in the rice crop classification.The grassland and the other crops that are spectrally similar to rice often create confusion for rice discrimination.Moreover, among the rice fields, variability exists due to the differences in crop planting time, local weather conditions and other factors [24].In this study, classification was performed on spectral data with and without using the phenology.The main differences in the results of the two approaches were the misclassification of grasslands and crop types.Using only the spectral features, the grasslands and other crops were misclassified as rice.The reason for the misclassification was that both the grassland and other crops had spectral properties similar to those of rice.On the other hand, the use of phenology along with the spectral features improved the classification significantly.The main reason behind the improvement was that the different crop types had different seasonal behavior, which is captured through the phenological variables.Additionally, the five phenological variables (base value, largest value, seasonal amplitude, large and small integral) are actually the statistical values of NDVI, meaning that they are calculated either by integrating, averaging, differencing or taking the maximum value of annual NDVI.These values are not dependent on the crop planting and harvesting date, thereby improving the classification by better capturing the crop growth profiles.
Discriminating a specific crop from the various vegetation types using a single image is a challenge [12,37].However, the images at the key crop growth stages effectively improve the classification accuracy.In the study area, the agricultural fields are generally small, fragmented and irregular.Therefore, a stepwise method has been established to identify the rice fields in a heterogeneous landscape using an object-based image analysis method.
This study demonstrates the applicability of multi-source remote sensing data (HJ-1A/B and MODIS) to accurate mapping of rice.The rice usually grows during the rainy season in a humid tropical climate where obtaining multi-temporal cloud-free optical images for the whole growing season is a constraint.This constraint becomes particularly critical when considering a large area.Therefore, optimal use of available cloud-free images is of great importance.This study demonstrates how the images from different sensors at peak crop growing stages help with achieving accurate rice classification.Additionally, the use of multi-source images increases the chance of getting more cloud-free images for the classification.
The top-down hierarchical classification approach yielded several advantages to discriminate rice fields.The general classes (e.g., water, forest, built-up, etc.) were classified by applying thresholds and then removed from further processing.The threshold condition for any class was adjusted to achieve higher classification accuracy without affecting the condition for the other classes.Along with the threshold condition, the nearest neighbor classification has also been applied for the initial classification of rice.The nearest neighbor classification considered the statistical distribution of the classes.It helped to identify the classes that could not be classified by the threshold condition.Moreover, the applied top-down hierarchical approach was logical to discriminate the rice class.
Use of object-based image analysis showed an advantage in identifying fragmented rice fields.Fragmentation of cropland is common in India [29].During the period between 2001 and 2010, the operational field size has reduced from 1.33 Ha down to 1.15 Ha [58].Land is becoming fragmented due to the absence of land use planning, rapid population growth, economic development, urbanization and the limited availability of arable land [13].This study demonstrated that the fragmented rice fields were well identified by the use of high resolution (30 m) HJ-1A/B images under the object-based image analysis framework.
In the recent study of Tian et al. [59] and Jia et al. [60,61], the number of base images used to produce the synthetic time series was much lower than the number used in this study, with six images in 12 years and one image in one year, respectively.In this study, three HJ CCD images of key rice growth stages were used as a base image for the ESTARFM fusion.When using different base images in the fusion process, the reported difference is minor between the synthetic and the actual images on the same date [59,62].Therefore, in spite of the three base images in our study, the accuracy of the blended time series was not significantly decreased.The base images were acquired at the beginning, middle and the end of the growing seasons, which helped the ESTARFM to capture the reflectance changes caused by phenology, and the synthetic time series reflected the actual changes in the NDVI trend [30].Additionally, the index-then-blend approach was adopted, which produces more accurate synthetic images [46].The adopted fusion technique yielded accurate synthetic NDVI.However, it is advised to use dense time series base images in order to achieve high fidelity in the synthetic images generated by ESTARFM.
Earlier works on rice mapping methods using optical images can be categorized into three groups.The first group considers individual cloud-free images and uses image statistics approaches, e.g., the unsupervised classifiers, like the self-organizing data analysis technique (ISODATA) [63,64], and supervised classifiers, like maximum likelihood [65] and support vector machine [66].
The second group considers time series images and uses different algorithms, like threshold-based and spectral matching techniques (SMT) [7,67].The third group is phenology-and pixel-based paddy rice mapping (PPPM) [15,17], where rice is identified for individual pixels based on the flooding signals of the rice transplanting phase by evaluating the differences between the Enhanced Vegetation Index (EVI)/NDVI and Land Surface Water Index (LSWI).The first two groups often generate maps that are difficult to compare for different regions, working groups and years, primarily due to the spectral heterogeneity, training sample selection, post-classification processing and the capabilities of the image interpreter [68].In comparison, the PPPM methods of the third category are less affected by these issues.The phenology-based methods can capture the growing stages of different crops and identify the unique signals according to crop calendar and management activities [34].In this study, rice growth stages were extracted by the analysis of 16-day MODIS NDVI and HJ CCD datasets, and they were used to improve the spectral-based mapping of rice.This study provides a rice mapping method combining phenology and multispectral satellite data (Table 6).To apply the proposed approach to other climatic regions and large regions, a few points need to be considered carefully, such as the appropriate selection of threshold conditions; parallel processing or cloud-computing is suggested for effective handling of large datasets [71]; and if possible, more high resolution images should be used.Increasing availability of medium and high spatial resolution satellite sensors like SPOT 6, Landsat 8, Sentinel-2 and G-F1 is increasing the observation frequency.An increase in good quality observations would increase the efficiency of the methodology, especially in the tropical regions, where cloud cover is common [72].
The identification of rice fields was affected by several potential factors.The problem of mixed pixels remains an issue.Some rice patches were too small to be identified with 30-m spatial resolution HJ-1A/B data.The widely-used, subtraction-based method could not be applied due to the unavailability of the SWIR (shortwave infrared) band in the HJ CCD sensor.The limited availability of cloud-free images for the whole growing season affected the rice identification.The vegetation dynamics could be better captured by using additional HJ-1A/B images during the sowing and harvesting period of rice.However, the use of phenology significantly reduces the necessity of the hyper-temporal image requirement for rice classification.

Conclusions and Future Works
This study presents an object-based image analysis approach for mapping of rice fields by the combined use of spectral and phenological features.The rice crop extent can provide important information for decision makers and government for proper management and monitoring of rice.Using the segmentation algorithm, image objects were created, which were then classified according to spectral, vegetation indices and phenological properties.It has been found that the proposed object-based image analysis approach yielded accurate classification results with an overall accuracy of 93% based on the validation samples.However, total rice areas were underestimated by ´34.53% when compared to the agricultural statistics.This study shows the importance of phenology for specific crop classification.This study also shows the applicability of just two season's multi-temporal high resolution images for rice crop classification.This is important because of the limited availability of cloud-free images.The approach presented can also be used to classify other crop types in different regions by making minor changes in classification threshold conditions and training samples.
Future work could be conduct on: (1) investigating the efficiency of phenological variables for fractional rice mapping in order to address the sub-pixel heterogeneity; (2) investigating the effects of fused time series NDVI on OBIA-based rice classification accuracy; (3) exploring the effects of image segmentation results on the OBIA-based rice classification; and (4) investigating more object features, such as textural and geometric features, for rice classification.

Figure 2 .
Figure 2. Calendar of the major crops in the study area.

Figure 2 .
Figure 2. Calendar of the major crops in the study area.

Figure 2 .
Figure 2. Calendar of the major crops in the study area.

Figure 3 .
Figure 3. Flowchart of object-based rice classification.ESTARFM, enhanced spatial and temporal adaptive reflectance fusion.

Figure 4 .
Figure 4. Raw and filtered NDVI time series for an agriculture pixel in the study region.

Figure 3 .
Figure 3. Flowchart of object-based rice classification.ESTARFM, enhanced spatial and temporal adaptive reflectance fusion.

Figure 3 .
Figure 3. Flowchart of object-based rice classification.ESTARFM, enhanced spatial and temporal adaptive reflectance fusion.

Figure 4 .
Figure 4. Raw and filtered NDVI time series for an agriculture pixel in the study region.Figure 4. Raw and filtered NDVI time series for an agriculture pixel in the study region.

Figure 4 .
Figure 4. Raw and filtered NDVI time series for an agriculture pixel in the study region.Figure 4. Raw and filtered NDVI time series for an agriculture pixel in the study region.

Figure 5 .
Figure 5. (a) HJ CCD NDVI of entire study area; (b) input MODIS NDVI; (c) synthetic MODIS NDVI; (d) HJ CCD NDVI nearest to the MODIS synthetic date.The red box in (a) indicates the enlarged area in (b), (c) and (d).

Figure 5 .
Figure 5. (a) HJ CCD NDVI of entire study area; (b) input MODIS NDVI; (c) synthetic MODIS NDVI; (d) HJ CCD NDVI nearest to the MODIS synthetic date.The red box in (a) indicates the enlarged area in (b-d).
of the season (SOS) Time for which the left edge has increased to 20% of the seasonal amplitude measured from the left minimum level Time of start of season End of the season (EOS) Time for which the right edge has decreased to 20% of the seasonal amplitude measured from the right minimum level Time of end of season Length of the season (LOS) Time from the start to the end of the season Total length of season Base value (BV) The average of the left and right minimum values Level of biomass Largest value (LV) Largest value between the start and end of the season Level of biomass Seasonal amplitude (SA) Difference between the maximum value and the base level Level of biomass Large integral (LI) Sum of all values from the season start to the end Net primary production Small integral (SI) Integral of the difference between the function defining season and the base value Net primary production 4.5.Classification -forest by the mean value of the red band.Again, the non-forest was classified into cropland and grassland by the mean NDVI value.The grassland tends to have high NDVI all year with less seasonal changes.On the other hand, croplands have variable NDVI with high values only in the peak growing seasons.Further, the NN supervised classification was applied to sub-divide the cropland class into rice and non-rice (grassland, built-up and other crops) classes.The mean value of the red, green and NIR bands, as well as phenological variables were utilized to perform the NN classification.Finally, the classes other than rice were re-assigned to one "non-rice" class in order to obtain a final rice and non-rice classification.To analyze the effect of the phenology on rice classification, two classification approaches were used: (1) using only the multi-spectral bands;

Figure 7 .
Figure 7. Classification scheme used for rice crop extraction.Doted boxes represent intermediate classes.

Figure 7 .
Figure 7. Classification scheme used for rice crop extraction.Doted boxes represent intermediate classes.
Remote Sens. 2016, 8, 479 10 of 20 phenological patterns.The concept of phenology is pointless in non-vegetated areas, which are also characterized by low NDVI variability over time.The hills and evergreen vegetation area show a long growing season (~176 days) with an earlier start and a later end.The variable used as a proxy for the green biomass, such as base value (BV) and the largest NDVI value (LV), shows lower values in the cultivated land, whereas higher values area recorded in the high biomass regions.The large integral (LI) value, which is related to net primary production[53], shows the highest values in the cropped land and forested area as compared to other land cover classes.The seasonal amplitude (SA) and small integral (SI), which is linked to plant biomass[57], varied according to the forest cover type, as evident from the southern part of the region.All of the phenological variables demonstrated typical patterns for the various vegetation classes and showed potential for crop identification.However, further research is required to fully utilize all of the variables for crop classification.

Figure 9 .
Figure 9. Classification results for different data combinations.(a) Classification results for October with spectral and spectral + phenological variable; (b) classification results for December with spectral and spectral + phenological variable; (c) classification results for March with spectral and spectral + phenological variable; (d) classification results for combined October, December and March with spectral and spectral + phenological variable.

Figure 9 .
Figure 9. Classification results for different data combinations.(a) Classification results for October with spectral and spectral + phenological variable; (b) classification results for December with spectral and spectral + phenological variable; (c) classification results for March with spectral and spectral + phenological variable; (d) classification results for combined October, December and March with spectral and spectral + phenological variable.

Figure 10 .
Figure 10.(a) Classification result differences between the spectral and spectral + phenological variable with respect to the October classification in Figure 9a; (b) classification result differences between the spectral and spectral + phenological variable with respect to the December classification in Figure 9b; (c) Classification result differences between the spectral and spectral + phenological variable with respect to the March classification in Figure 9c; (d) classification result differences between the spectral and spectral + phenological variable with respect to the combined October, December and March classification in Figure 9d.

Figure 10 .
Figure 10.(a) Classification result differences between the spectral and spectral + phenological variable with respect to the October classification in Figure 9a; (b) classification result differences between the spectral and spectral + phenological variable with respect to the December classification in Figure 9b; (c) Classification result differences between the spectral and spectral + phenological variable with respect to the March classification in Figure 9c; (d) classification result differences between the spectral and spectral + phenological variable with respect to the combined October, December and March classification in Figure 9d.

Figure 11 .
Figure 11.Spatial distribution of the derived rice cultivated area.

Figure 11 .
Figure 11.Spatial distribution of the derived rice cultivated area.

Table 1 .
Satellite datasets used for the study.

Table 5 .
Difference between estimated rice area and statistics (thousand ha).

Table 6 .
Comparison of some of the existing phenology-based paddy rice mapping algorithms.