NESEA-Rice10: high-resolution annual paddy rice maps for Northeast and Southeast Asia from 2017 to 2019

An accurate paddy rice map is crucial for ensuring food security, particularly for Southeast and Northeast Asia. MODIS satellite data are useful for mapping paddy rice at continental scales but have a mixed-pixel problem caused by the coarse spatial resolution. To reduce the mixed pixels, we designed a rule-based method for mapping paddy rice by integrating time series Sentinel-1 and MODIS data. We demonstrated the method by generating annual paddy rice maps for Southeast and Northeast Asia in 2017–2019 (NESEA-Rice10). We compared the resultant paddy rice maps with available agricultural statistics at subnational levels and existing rice maps for some countries. The results demonstrated that the linear coefficient of determination (R²) between our paddy rice maps and agricultural statistics ranged from 0.80 to 0.97. The paddy rice planting areas in 2017 were spatially consistent with the existing maps in Vietnam (R²=0.93) and Northeast China (R²=0.99). The spatial distribution of the 2017–2019 composite paddy rice map was consistent with that of the rice map from the International Rice Research Institute. The paddy rice planting area may have been underestimated in the region in which the flooding signal was not strong. The dataset is useful for water resource management, rice growth, and yield monitoring. The full product is publicly available at https://doi.org/10.5281/zenodo.5645344 (Han et al., 2021a). Small examples can be found from the following DOI: https://doi.org/10.17632/cnc3tkbwcm.1 (Han et al., 2021b).

1 Introduction

Rice is one of the world's main food sources, accounting for approximately 12 % of the global cropland area (Zhang et al., 2018; Singha et al., 2019). Approximately 90 % of the world's rice is produced in Asian countries (Chen et al., 2012; Yeom et al., 2021). Rice provides food for over 50 % of the world's population (Minasny et al., 2019). The consumption of rice increases as the world's population increases. Additionally, approximately $\mathrm{1}/\mathrm{10}$ of CH₄ emissions in the atmosphere come from methane emissions from rice paddies (Ehhalt et al., 2001; Xin et al., 2017; Zhang et al., 2020). Rice agriculture is significant in food security, water resource security, disease transmission, and environmental sustainability (Clauss et al., 2018b; Li et al., 2020; Park et al., 2018). An accurate planting area and spatial distribution information are the basis for monitoring paddy rice growth and predicting yield. However, few spatial maps of paddy fields at continental scales exist (Li et al., 2020; Singha et al., 2019). Therefore, it is necessary to produce a paddy rice map dataset with high spatial resolution.

Many methods for mapping rice have been developed based on different remote sensing data, including (1) machine learning classifiers (e.g., random forest and support vector machines), (2) phenology-based classifiers, (3) rule-based algorithms, and (4) the time series algorithm approach (Dong et al., 2016a, b; Bazzi et al., 2019; Dong and Xiao, 2016; Luo et al., 2020b, a; Nelson et al., 2014; Phung et al., 2020; Minasny et al., 2019; Shew and Ghosh, 2019; Xiao et al., 2006; Zhan et al., 2021). Satellite image sources include MODIS, Landsat, Sentinel, RADARSAT, and PALSAR (Dong and Xiao, 2016; Shao et al., 2001; Singha et al., 2019; Zhou et al., 2016). Many studies have demonstrated that phenology-based classifiers using MODIS data are useful for mapping paddy rice at continental scales (Dong et al., 2016b; Xiao et al., 2006; Zhang et al., 2020). The transplanting period of rice is a distinct characteristic used for distinguishing rice from other crops or land use types. For example, Xiao et al. (2006) mapped paddy rice at continental scales (South Asia and Southeast Asia, SE Asia) using the phenological characteristics in the period of flooding and transplanting. Additionally, this method was successfully applied in other large regions (Xin et al., 2020; Zhang et al., 2017, 2020). The International Rice Research Institute (IRRI) extracted the distribution of paddy rice for Asia (Nelson and Gumma, 2015). However, the paddy rice maps generated using MODIS data contain a large number of mixed pixels caused by the coarse spatial resolution (500 m) (Dong et al., 2015, 2016b; Shew and Ghosh, 2019), particularly in hilly areas (Z. Liu et al., 2019). The mixed land cover types within MODIS pixels can affect the accuracy of the rice map (Sun et al., 2009). Fine spatial resolution images, including Landsat TM/ETM+/OLI, Huan Jing (HJ), and Sentinel-2 images, are also used for mapping paddy rice. Some previous studies have shown that rice maps generated from Landsat images have relatively high accuracy (Dong et al., 2016a; Torbick et al., 2017). However, they are only suitable for relatively small study areas in which cloud cover is minimal and not for continental scales (Ramadhani et al., 2020; Torbick et al., 2017). In contrast to optical satellite images, synthetic aperture radar (SAR) data are unaffected by clouds (Park et al., 2018). Moreover, SAR data have the special characteristics of backscatter changes during the growth of paddy rice (Bazzi et al., 2019; Clauss et al., 2018b; Nguyen et al., 2016; Phung et al., 2020; Planque et al., 2021). For example, Singha et al. (2019) mapped the rice map for Bangladesh based on Sentinel-1 data and random forest classifiers with good accuracy. Although paddy rice has been mapped in several studies using SAR data, they are still difficult to use widely over large areas because of the lack of a large number of ground truth samples (Clauss et al., 2018b; Minh et al., 2019; Nguyen et al., 2016; Phung et al., 2020; Minasny et al., 2019; Zhan et al., 2021; Zhang et al., 2018). Because the average area of crop fields in many regions in Asia is less than half a hectare (Maclean et al., 2013), it is critical to generate paddy rice maps with higher spatial resolution at continental scales than past efforts with MODIS.

Optical remote sensing images and SAR data have complementary information (Park et al., 2018; Wang et al., 2015). The combination of optical and SAR images can provide opportunities for mapping paddy rice with a few mixed pixels and a high spatial resolution at continental scales. MODIS data have the advantage of high temporal resolution, which reduces cloud problems and provides valuable spectral information for identifying paddy rice. Sentinel-1 SAR data with a high spatial resolution (10 m) provide backscatter information for different land types. Therefore, the integration of MODIS and SAR images may solve the mixed-pixel issue to a great degree and enable the production of more reliable paddy rice maps than those based only on MODIS images (Dong and Xiao, 2016; Park et al., 2018; Torbick et al., 2010; Wang et al., 2015). We take advantage of both MODIS and SAR strengths to map paddy rice fields at a large scale.

Thus, we aim to improve the MODIS-based method for mapping paddy rice fields by integrating Sentinel-1 SAR data to reduce mixed-pixel effects. Then we use the method to generate paddy rice maps in 2017–2019 for SE Asia and Northeast Asia (NE Asia). The map products will be useful for scientific communities and stakeholders for many purposes.

2 Materials

2.1 Study area

The study areas were NE and SE Asia. NE Asia is composed of Northeast China (Liaoning, Jilin, and Heilongjiang province), the Democratic People's Republic of Korea, the Republic of Korea, and Japan (Dong et al., 2016b; Yeom et al., 2021). The main paddy-rice-producing regions in NE Asia are concentrated in the plain in Northeast China, the western plain of the Korean Peninsula, and the alluvial plains around the Japanese islands. In SE Asia, the countries where rice is planted intensively include Indonesia, Thailand, Vietnam, Myanmar, the Philippines, Malaysia, and Myanmar. SE Asia cultivates approximately 30 % of the world's rice (Bridhikitti and Overcamp, 2012; Huke and Huke, 1997). The dense planting areas of rice in SE Asia are located in valleys and deltas, such as the Red River delta in northern Vietnam and the Mekong River delta in southern Vietnam (Clauss et al., 2018a; Phung et al., 2020). The Mekong delta produces more than half the rice in Vietnam (Bouvet and Le Toan, 2011). The main rice cropping system in NE Asia is single rice (Dong et al., 2016b). By contrast, three rice cropping systems are dominant in SE Asia: single rice, double rice, and triple rice (Laborte et al., 2017). Because climate and crop calendars vary across SE and NE Asia, the study area was classified into eight refined agroecological zones based on temperature, seasonal precipitation, and farming practices from previous studies (Oliphant et al., 2019; Suepa et al., 2016). The zones were further subdivided into 41 regions for classification (Fig. 1).

Figure 1 Agroecological zones and 100 m radius sample blocks in SE and NE Asia.

2.2 Data source

2.2.1 Satellite imagery

We acquired the time series MOD09A1 images from the Google Earth Engine (GEE) data catalog (https://developers.google.com/earth-engine/datasets/, last access: 18 June 2021). The 8 d composite MOD09A1 provides seven surface reflectance bands (red, blue, green, near-infrared, and shortwave infrared 1–2 bands) at 500 m spatial resolution. We used the blue band (≥0.2) to remove cloudy pixels (Chen et al., 2012; Xiao et al., 2006). We projected MODIS data using the WGS1984 coordinate system. Additionally, we collected SAR data with 10 m resolution in interferometric wide-swath mode from the Sentinel-1 satellite (Torres et al., 2012). To maximize the frequency of observations, we used the Level-1 ground range detected (GRD) product in ascending and descending orbits. We pre-processed the GRD data with VH (Vertical transmit Horizontal receive) and VV (Vertical transmit Vertical receive) polarization that we acquired (e.g., calibration and geocoding) using the Sentinel-1 Toolbox on the GEE (https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD, last access: 18 June 2021) (Singha et al., 2019; C. Liu et al., 2020). Following this, we applied a filter with a boxcar kernel (30×30) moving to reduce speckle noise (Minasny et al., 2019). However, the rugged terrain and side-looking SAR imaging geometry caused radiometric distortion (Mullissa et al., 2021; Vollrath et al., 2020). To increase the quality of observations, we used the physical reference model (volume) proposed by Vollrath et al. (2020) to make a radiometric slope correction for Sentinel-1 on the GEE. Figure S1 in the Supplement shows an example of the effect of slope correction in SE Asia using the model. The code can be found at https://github.com/ESA-PhiLab/radiometric-slope-correction (last access: 18 June 2021).

2.2.2 Terrain data

We generated digital elevation model (DEM) data from the Shuttle Radar Topography Mission (SRTM) Version 4 (Reuter et al., 2007). The spatial resolution of the DEM was 90 m × 90 m. We acquired the DEM data and calculated the slope map from the DEM on the GEE platform (Table 1).

Table 1 Detailed information about the data used in this study.

2.2.3 Forest land

We extracted the forest land mask from the Global PALSAR Forest Map in 2017 (Table 1). The Global PALSAR Forest Map (25 m spatial resolution) was generated by the Japan Aerospace Exploration Agency (JAXA) (Shimada et al., 2014). Pixels with a forest area larger than 0.5 ha and forests covering over 10 % of the pixel area were defined as forest pixels (Shimada et al., 2014).

2.2.4 Wetland

We extracted the distribution of wetland from the GlobeLand30 dataset in 2020. GlobeLand30 is available from the National Geomatics Center of China (Table 1). This product at 30 m spatial resolution with high accuracy was generated using Landsat, Chinese HJ-1, and GF-1 satellite images (http://www.globallandcover.com, last access: 18 June 2021) (Chen et al., 2015).

Finally, we resampled all the raster data to 10 m to match the spatial resolution of Sentinel-1.

2.2.5 Agricultural statistics

We collected annual rice planting area census data at the subnational level (state, province, city, prefecture, or county) from the available statistical yearbooks of various countries. The agricultural statistics were provided by agricultural statistical offices. The areas in the statistics data were converted into hectares (ha). Detailed information about the collected agricultural statistics in this study is presented in Table 1.

2.2.6 Existing rice maps

We collected the existing publicly available rice maps from three sources: (1) the 500 m spatial resolution paddy rice map with high accuracy in Southern China in 2017 that was generated using the phenology- and pixel-based algorithm from MODIS data (Xin et al., 2020), (2) the High-Resolution Land Use and Land Cover (HRLULC) map for Vietnam in 2017 (Hashimoto et al., 2014) with 10 m spatial resolution generated using multiple remote sensing data sources, and (3) the 500 m resolution rice maps of Asia obtained from the IRRI (Nelson and Gumma, 2015), which were mainly derived from MODIS data. We compared these existing products with our paddy rice maps.

2.3 Methodology

2.3.1 Analyzing the characteristics of spectral indices from paddy rice

There are three growing stages for paddy rice: transplanting, growing, and post-harvest periods (Singha et al., 2019). Flooding signals in the transplanting period are unique characteristics that distinguish paddy rice from other crops (Clauss et al., 2016; Dong et al., 2016b; Sun et al., 2009). The color combination of MODIS images (R/G/B = band7/band2/band1) in the transplanting stage of the paddy rice field has a prominent tone (Fig. S5). We calculated the land surface water index (LSWI) and enhanced vegetation index (EVI) for each image:

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{LSWI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{SWIR}}}{{\mathit{\rho }}_{\mathrm{NIR}}+{\mathit{\rho }}_{\mathrm{SWIR}}}},\text{(2)}& {\displaystyle }\text{EVI}=\mathrm{2.5}\times {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{RED}}}{{\mathit{\rho }}_{\mathrm{NIR}}+\mathrm{6}\times {\mathit{\rho }}_{\mathrm{RED}}-\mathrm{7.5}\times {\mathit{\rho }}_{\mathrm{BLUE}}+\mathrm{1}}},\end{array}$$

where ρ_SWIR, ρ_NIR, ρ_RED, and ρ_BLUE are values of band 6, band 2, band 1, and band 3, respectively. Note that we chose EVI instead of the normalized difference vegetation index (NDVI) for paddy rice identification because NDVI is more sensitive to atmospheric contamination and has a saturation issue (Zhang et al., 2015). Figure 2 shows the standard temporal profile of the EVI and LSWI of different paddy rice planting systems (single, double, and triple) at three typical sites. When LSWI plus 0.05 is larger than EVI, it indicates that the paddy rice is in the transplanting period. Both the EVI and LSWI values increase after paddy rice is transplanted. EVI values are higher than LSWI values because the fields are fully covered by the rice canopy. The EVI decreases during the post-harvest period. Both double paddy rice and triple paddy rice have flooding signals. The above phenomena are consistent with previous studies (Dong et al., 2016b; Minh et al., 2019; Shew and Ghosh, 2019; Xiao et al., 2006). The profiles of the LSWI and EVI of some land cover types (e.g., water, urban, and forest) are different from paddy rice (Fig. S3). Therefore, color gradations and the relationship between EVI and LSWI are useful for extracting phenological information of paddy rice. In this study, we did not smooth the EVI and LSWI time series. LSWI changes under different dry and wet conditions and smoothing the EVI and LSWI time series may eliminate the true paddy rice flood signal (L. Liu et al., 2020). Therefore, we did not reconstruct the EVI and LSWI datasets.

Figure 2 Temporal profile analysis of EVI, LSWI, VV, VH, and VH/VV ratio from three typical paddy rice sites with different latitudes during 2017–2020: (a) single paddy rice in Northeast China (46.804562^∘ N, 131.896961^∘ E), (b) double paddy rice in the Philippines (15.630637^∘ N, 120.864242^∘ E), and (c) triple paddy rice in Indonesia (8.6211997^∘ S, 116.141157^∘ E). The shaded areas indicate the standard deviation.

2.3.2 Analyzing the characteristics of backscatter coefficients from paddy rice

The backscatter coefficients change as paddy rice grows and develops. Paddy rice fields appear as a black area in the VH image on the transplanting date (Fig. S2) because the water (flood) in the transplanting period decreases the VH backscatter coefficient values (Dineshkumar et al., 2019; Torbick et al., 2017). The VH and VV backscatter coefficients have a local minimum value during the transplanting period in all reference paddy rice fields (Fig. 2). After transplanting, the VH backscatter coefficients increase as the paddy rice grows and reaches a peak at the heading stage (Zhan et al., 2021; Zhang et al., 2018). The VH backscatter coefficients decrease after the rice harvest stage (Phung et al., 2020; Singha et al., 2019; Torbick et al., 2017). Additionally, paddy rice has consistent temporal behavior in the VH/VV ratio and VH. The profiles of the dynamic backscatter coefficients of some land cover types (e.g., water, urban, and forest) are different from those of paddy rice (Fig. S4). Therefore, color gradations and the time series of backscatter coefficients are useful for identifying paddy rice phenology information (Yulianto et al., 2019; Phung et al., 2020; Zhan et al., 2021).

2.3.3 Sample blocks collected for extracting phenological parameters

Paddy rice in SE and NE Asia is cultivated using diverse cropping systems because of the climate and other natural conditions (Dong et al., 2016a; Laborte et al., 2017; Nelson et al., 2014; Shew and Ghosh, 2019). With reference to previous studies (Clauss et al., 2016; Gumma et al., 2014; C. Liu et al., 2020; Phan et al., 2019; Phung et al., 2020), we acquired information about the flooding signal period and length of the growing season for each subzone using sampling-based information. We selected sample blocks that were distributed over the different rice-growing zones across SE and NE Asia. Each block was a polygon with a radius of 100 m. We collected the sample blocks according to multiple rules (Clauss et al., 2016; Dong et al., 2016a; Fikriyah et al., 2019; Singha et al., 2019). First, the time series of the backscatter coefficients and vegetation index of the mean values from all pixels in each sample block were consistent with the phenological characteristics of paddy rice (Sect. 2.3.1 and 2.3.2). Second, the sample blocks were also digitized using Google Earth or Sentinel-1/2 images using visual interpretation referring to previous studies (Dong et al., 2016a; Zhang et al., 2015). Third, we also used existing rice maps and calendar information as complementary information (Laborte et al., 2017; Maclean et al., 2013). Note that not all Google Earth or Sentinel-2 images were available throughout SE and NE Asia. We collected a total of 438 sample blocks and 504 sample blocks using the above rules for SE and NE Asia, respectively (Fig. 1). These blocks covered most paddy rice fields in the study areas. We generated mean backscatter coefficients and vegetation index time series profiles for each block. Following this, we manually extracted the paddy rice growth and phenological parameters based on the backscatter time series characteristics. Finally, we obtained the phenological information for each subzone from the sample blocks (Gumma et al., 2014; C. Liu et al., 2020). Although there may be some limitations in extracting phenological parameters for zones using random samples, it may be one of the most effective approaches currently available (Clauss et al., 2016; Gumma et al., 2014; Han et al., 2021c; Li et al., 2020; Phan et al., 2019; Phung et al., 2020).

Figure 3 Flowchart for mapping paddy rice in SE and NE Asia using multiple data.

2.3.4 Algorithm for identifying paddy rice fields

We used a rule-based method to map paddy rice and produce annual paddy rice maps for SE and NE Asia in 2017–2019 at 10 m resolution (NESEA-Rice10) using the phenological features of paddy rice (Fig. 3). The steps for generating the paddy rice maps are as follows.

1. Detect the flooding area of paddy rice. The key features used to identify paddy rice are the flooding signals in the transplanting phase (Dong and Xiao, 2016). We used LSWI $+\mathrm{0.05}>$ EVI and LSWI <0.45 to extract the flooding signals. This rule has been used to successfully map potential paddy rice fields over large areas (Sakamoto et al., 2009; Xiao et al., 2006; Zhou et al., 2016). Because paddy rice flooding signals occur over a short period, non-permanent flooding (e.g., persistent water bodies and fishponds) should be removed (Nelson et al., 2014; Zhang et al., 2015). We removed pixels that had more than 20 composite periods identified as flooding signals during a year.

In addition to the optical MODIS-based LSWI–EVI relationship approach, we also applied the minimum value of VH data in the transplanting stage to identify flooding signals, as suggested in previous studies (Clauss et al., 2018b). VH has a higher sensitivity in paddy rice growth stages than VV polarization (Inoue et al., 2020; Nguyen et al., 2016; Wakabayashi et al., 2019). However, the minimum value of VH in different regions is different because Sentinel-1 data are affected by the incidence angle (ranging from approximately 30 to 45^∘) (Fig. S6) (Phung et al., 2020; Singha et al., 2019; Zhang et al., 2018). Currently, it is still challenging to normalize the incidence angle over a large area on the GEE. The VH value of the water surface changes continuously from −21 to −34 dB as the incidence angle changes (Phung et al., 2020). To reduce the incidence angle effect on the Sentinel-1 images, we considered −20 dB as the conservative baseline threshold to achieve the minimum VH value. In previous studies, researchers also proved the effectiveness of this threshold for identifying flooding signals in paddy fields (Clauss et al., 2018b; Nguyen et al., 2015, 2016; Nguyen and Wagner, 2017; Zhang et al., 2018). To further improve the accuracy of flooding signal extraction, we fine-tuned the baseline threshold of VH for each subzone based on the histogram of sample blocks collected in Sect. 2.3.3. We considered the pixels that met all the above conditions as flooding signals.

2. The EVI of paddy rice increased rapidly after the transplanting period because of the increasing numbers of leaves and biomass (Chen et al., 2012). Therefore, we removed pixels with a maximum EVI value of less than 0.4, as suggested in previous studies (Kontgis et al., 2015; Sakamoto et al., 2009).

3. Moreover, the coefficient of variation (CV) has been proven to be an effective indicator for distinguishing crop types and non-cropland (Huang et al., 2021; C. Liu et al., 2020; Rose et al., 2021; Whelen and Siqueira, 2018). The VH backscatter coefficient of crops, particularly of paddy rice, has a larger time series variation range than non-agricultural land (e.g., urban and water) (Fig. S4). Based on this time series characteristic, paddy rice may be identified from different land surfaces. Therefore, we removed pixels with CV values greater than −0.4 and less than −0.1 calculated using VH during the growth of paddy rice (Rose et al., 2021; Whelen and Siqueira, 2018). The threshold was determined from the histogram of the paddy rice sample blocks in the study area. The histograms can be found in Fig. S7. We measured the CV_VH using the temporal mean (mean) and standard deviation (SD) of the time series of VH during the paddy rice growth period:

$$\begin{array}{}\text{(3)}& {\text{CV}}_{\mathrm{VH}}={\displaystyle \frac{\text{SD}}{\text{mean}}}.\end{array}$$

4. We used the mask with slopes larger than 5^∘ to remove steep terrain; it is unsuitable to plant paddy rice on sloping land (Sun et al., 2009).

5. We used the PALSAR-based forest map in 2017 as a mask (Wang et al., 2015; Zhang et al., 2017).

6. We used the water mask from GlobalLand30 in 2020 referring to the study of Zhang et al. (2015) to reduce the misclassification of paddy rice as it is challenging to extract paddy rice from wetlands because paddy rice and wetland have similar characteristics to flooding signals (Zhang et al., 2015; Zhou et al., 2016). 7. We classified pixels that met all of the above rules as paddy rice. Following this, we deleted small isolated pixels (connected components less than 12 pixels) to remove the “salt and pepper” effect in the classification (https://catalog.data.gov/dataset/global-food-security-support-analysis-data-gfsad-cropland-extent-2015-southeast-and-northe, last access: 18 June 2021).

The single cropping system for paddy rice identification is not ideal because of the difference in paddy rice cultivation time in some regions of SE Asia (Fikriyah et al., 2019; Shew and Ghosh, 2019). Therefore, we combined all paddy rice fields identified at different times into the annual map. We applied the improved method to generate the annual paddy rice maps for SE and NE Asia in 2017–2019. Please note that the method we improved may not extract rice fields (e.g., rain-fed paddy rice and upland rice) if flooding signals are not available (Xiao et al., 2006; Zhang et al., 2017).

2.4 Accuracy assessment

It is challenging to evaluate the accuracy of the classification at continental scales (Xiao et al., 2006; Zhang et al., 2020). We used two strategies to evaluate the paddy rice maps as accurately as possible. First, we compared the available agricultural statistics on a subnational level in some countries (Table 1). Referring to the study of Xiao et al. (2006), we calculated the annual area of paddy rice based on paddy intensity. The paddy intensities of countries in NE Asia, Myanmar, Vietnam, and the Philippines were 1, 1.4, 2.2, and 2, respectively. Second, we compared the spatial consistency between our classification results and existing rice maps (Table 1). We used the coefficient of determination (R²) to measure the consistency between our paddy rice maps, agricultural statistics, and existing products.

$$\begin{array}{}\text{(4)}& R^{\mathrm{2}}={\displaystyle \frac{{\left({\sum }_{i=\mathrm{1}}^n\left(x_i-{\stackrel{\mathrm{‾}}{x}}_i\right)\times \left(k_i-{\stackrel{\mathrm{‾}}{k}}_i\right)\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(x_i-{\stackrel{\mathrm{‾}}{x}}_i\right)}^{\mathrm{2}}\times {\sum }_{i=\mathrm{1}}^n{\left(k_i-{\stackrel{\mathrm{‾}}{k}}_i\right)}^{\mathrm{2}}}},\end{array}$$

where n is the total number of administrative units, x_i represents the mapped paddy rice areas, $\stackrel{\mathrm{‾}}{x_i}$ is the corresponding mean value, k_i represents the agricultural statistics or areas from existing rice maps, and $\stackrel{\mathrm{‾}}{k_i}$ is the corresponding mean value.

Figure 4 Comparison of the resultant annual paddy rice areas and the agricultural statistics at the subnational level in different countries from 2017 to 2019. The marginal kernel density plot above or to the right of each scatterplot shows the distribution of the data in one dimension.

3 Results

3.1 Comparison of the classification with agricultural statistics

The paddy rice maps in SE and NE Asia in 2017–2019 are presented in Figs. S8 and S9, respectively. We calculated the annual paddy rice area using the pixel number approach for each administrative unit. The estimated annual rice paddy areas were significantly correlated with the agricultural statistics at subnational levels. The resultant paddy rice maps and the agricultural statistics had relatively high correlations in Northeast China (R² ranged from 0.82 to 0.89, p<0.01) (Fig. 4a). The paddy areas in Changchun, Jilin, and Tonghua were underestimated. This is mainly because of the lack of available satellite data. When we excluded the three cities, R² ranged from 0.85 to 0.97, with significant correlations. Additionally, there were significant correlations between the paddy rice maps and agricultural statistics in the Republic of Korea (Fig. 4b), but the results underestimated the paddy rice area (R² ranged from 0.80 to 0.82, p<0.01). The main reason may be that many small rice fields were situated in narrow valleys in the mountains (Dong et al., 2016b; Peng et al., 2011). The correlations were high in most counties in Japan (R² ranged from 0.89 to 0.93, p<0.01). In Hokkaido, the paddy rice areas were overestimated (Fig. 4c), as they were by Zhang et al. (2018). The resultant paddy rice areas were consistent with the agricultural statistics in Myanmar, with R² ranging from 0.91 to 0.94 (Fig. 4d), and in Vietnam, with R² equal to 0.97 in the 3 years (Fig. 4e). R² between the paddy rice maps and the statistical data ranged from 0.81 to 0.87 in the Philippines, but rice areas in some provinces were underestimated (Fig. 4f). Cloud contamination may be a major reason for the underestimation (Peng et al., 2011; Xiao et al., 2006). The spatial distribution of paddy rice was visually consistent with that of the higher spatial resolution images in some typical testing regions (Fig. S10).

Figure 5 Comparison of the annual paddy rice area between our classification and existing datasets at the subnational level in Northeast China (a) and Vietnam (b). The marginal kernel density plot above or to the right of each scatterplot shows the distribution of the data in one dimension.

3.2 Comparison of the classification with other annually available rice maps

We further compared the resultant rice maps with existing rice maps at the subnational level. The annually available datasets included the MODIS-based rice paddy map with 500 m resolution for Northeast China in 2017 and the JAXA-derived rice map with a 10 m resolution for Vietnam in 2017 (Sect. 2.2.6). The paddy rice area statistics from our maps and existing products significantly correlated with R²=0.99 (p<0.01) for Northeast China and R²= 0.93 (p<0.01) for Vietnam (Fig. 5). We note that the paddy rice area in our maps was smaller than that in the MODIS-based product for Northeast China. The main reason may be that the 500 m resolution MODIS-based paddy rice map had a large number of mixed pixels. Although the spatial patterns of our maps were consistent with the MODIS-based products (Figs. S11 and S12), our maps contained more details with fewer mixed pixels (Fig. 9a–f). Additionally, the detailed information at the field scale was consistent for both our maps and the JAXA-derived paddy rice map (Fig. 9g–r). Overall, the comparison of the classification with existing products confirmed the reliability of the paddy rice maps that we generated.

Figure 6 Comparisons between the paddy rice area in our study and IRRI dataset in SE Asia at the (a) national and (b) county levels. The marginal kernel density plot above or to the right of each scatterplot shows the distribution of the data in one dimension.

3.3 Comparison of the composite paddy rice map and IRRI dataset

The composite paddy rice map is a mosaic of rice planting areas in 3 years (2017–2019) where rice has been detected in 1 or more years. We compared the composite paddy rice areas with IRRI products at the national and subnational levels in SE and NE Asia. The results demonstrated that the correlations between them were significant at both levels (R² ranged from 0.73 to 0.80) (Fig. 6). The paddy rice area based on the IRRI product was higher than our results. The main reason may be that our method reduced the mixed pixels in the paddy rice map and that the IRRI product from MODIS overestimated the area, as in previous studies (Fig. S13) (Chen et al., 2012; Li et al., 2020; Nelson and Gumma, 2015). The difference between the two datasets may also be partly caused by the inconsistent epoch composite years in the two datasets. Moreover, although the distribution of paddy rice was consistent between our results and the IRRI product, there were regional differences. Our paddy rice map underestimated the rice area in Thailand (Fig. 6a), which may be because the rice planted in eastern Thailand has no obvious strong flooding signals (Zhang et al., 2020). Despite the various spatial resolutions and different years in the rice paddy data, the intercomparison verified the accuracy of our dataset.

3.4 Spatial patterns of paddy rice areas

In NE Asia, paddy rice fields are primarily cultivated in the longitude range from 123 to 134^∘ E and latitude range from 45 to 48^∘ N (Fig. 7). In Northeast China, paddy rice is mainly cultivated in Heilongjiang (Sanjiang Plain) and Liaoning provinces. Paddy rice in the Democratic People's Republic of Korea and the Republic of Korea is distributed in the western coastal plains. Some paddy rice fields are located in narrow valleys. Rice fields in Japan are mainly on the coastal alluvial plains.

Figure 7 Spatial distribution of classified composite paddy rice with a 10 m spatial resolution in NE Asia during 2017–2019. The curves represent the relative change rate in the distribution of the number of paddy rice pixels along the longitude and latitude gradients.

Paddy rice is generally planted in the plains and deltas of rivers in SE Asia in the latitude range from 10 to 21^∘ N and the longitude range from 94 to 106^∘ E (Fig. 8). For example, the Mekong Basin and Hong River delta are typical main rice-growing areas in Vietnam. Paddy rice in the Philippines is mostly grown in the northern plains and scattered in the southern Philippines. Most of the rice cultivation in Malaysia is in the northwest corner of the peninsular region.

Figure 8 Spatial distribution of classified composite paddy rice with 10 m spatial resolution in SE Asia during 2017–2019. The curves represent the relative change rate in the distribution of the number of paddy rice pixels along the longitude and latitude gradients.

Figure 9 Visual comparison of our paddy rice maps and existing products in typical regions in 2017: (a–c, g–i, m–o) classification using our method. The example data for (a, g, m) can be found at the following link: https://doi.org/10.17632/cnc3tkbwcm.1 (Han et al., 2021b); see example 01–03. (d–f) MODIS-based paddy rice fields and (j–l, p–r) JAXA-derived rice fields are also shown.

4 Discussion

MODIS data were useful for mapping paddy rice at continental scales using combined EVI and LSWI analysis. Most paddy rice fields were fragmented in Asia (Li et al., 2020; Lowder et al., 2016). Therefore, it is difficult to solve the intra-class temporal variability of paddy rice pixels caused by the coarse resolution of 500 m (Dong et al., 2016b; Xiao et al., 2006). Mixed pixels may cause an overestimation of the rice cultivation areas (Nelson and Gumma, 2015). We improved the MODIS-based approach by incorporating Sentinel-1 data and used the approach to identify paddy rice fields in SE and NE Asia for 2017–2019. Reducing the mixed-pixel problem is the key point of the improved paddy rice mapping method. Compared with the paddy rice maps acquired from existing MODIS-based products, our classification provides more information about field details with a higher spatial resolution (10 m) (Fig. 9). Therefore, the integration of MODIS and Sentinel-1 data makes it possible to improve the accuracy of mapping paddy rice at continental scales.

Although our paddy rice maps are consistent with existing products, some uncertainty sources still affect the mapping results. First, identifying small paddy rice fields in hilly regions is challenging for MODIS data, which will lead to an underestimation of the area of paddy rice fields (Dong and Xiao, 2016; Zhang et al., 2015). For example, the rice planting area is smaller than the agricultural statistics in the mountainous provinces of the Republic of Korea (Fig. 10). The classification method relies on rice paddies containing irrigation water during transplanting stages. Therefore, rain-fed paddy rice and upland rice may not be detected because of the unavailability of flooding signals (Zhang et al., 2017). The main reason for the underestimation of the rice area in eastern Thailand may be that the flooding signal of rice was not detected, which has also been mentioned in previous studies (Bridhikitti and Overcamp, 2012; Guo et al., 2019; Zhang et al., 2020). Although MODIS data with a high temporal resolution was used in our method, the accuracy of rice maps is still affected by cloud contamination (Fig. 11a–b) (Dong and Xiao, 2016). Missing observations in Sentinel-1 data would lead to noteworthy omission errors (Fig. 11c–d). In addition, both the thresholds of different indicators and phenological information extracted by sample blocks may affect the accuracy of the classification (Dirgahayu and Parsa, 2019; Jeong et al., 2012; Li et al., 2020; Yeom et al., 2021). Finally, uncertainties in other land cover products used in this study may also affect the accuracy of the classification.

Figure 10 Estimated distribution of paddy rice in 2017 in mountainous regions in South Korea: (a) flooding signal based on MODIS indices, a (b) paddy rice map generated by our method, and (c) mountainous landscapes acquired from © Google Earth.. The example data for b can be found at the following link: https://doi.org/10.17632/cnc3tkbwcm.1 (Han et al., 2021b; see example 04).

Figure 11 Spatial distribution map of good-quality observation numbers during 2017 to 2019 for (a1–a3) MODIS images in Northeast Asia, (b1–b3) MODIS images in Southeast Asia, (c1–c3) Sentinel-1 images in Northeast Asia, and (d1–d3) Sentinel-1 images in Southeast Asia.

Under the combined effects of climate change and human activities, such as frequent extreme disasters, population growth, and urban expansion, knowing the spatial distribution of paddy rice is important for food security. The potential applications of the dataset include (1) improving paddy rice yield prediction accuracy, as crop masks are the basis for paddy rice yield prediction, and previous studies have demonstrated that the accuracy of crop masks affects the accuracy of yield prediction (J. Liu et al., 2019; Zhang et al., 2019); (2) assessing damage to agriculture from extreme hazards as floods are one of the major natural disasters in Southeast Asia, and high-resolution paddy rice maps will improve the accuracy of the area and yield loss estimates for flooded farmland (Phan et al., 2019); and (3) estimating greenhouse-relevant methane emissions. Paddy rice is an important source of methane in the atmosphere (Redeker et al., 2000). Accurate paddy rice maps and crop intensity maps facilitate the estimation of methane emissions (Zhang et al., 2020). In addition, paddy maps are helpful in making land-use decisions for the government.

Recently, as more Sentinel-2 images with higher resolutions have become available, combining Sentinel-2 and other satellite images have improved the temporal resolution of the data. For example, the Harmonized Landsat and Sentinel-2 (HLS) project provide images with 2–3 d at 30 m spatial resolution by combining Landsat 8 satellite and Sentinel-2 satellite data (https://hls.gsfc.nasa.gov/, last access: 5 November 2021). Zhang et al. (2021) mapped the global cropping intensity with a high spatial resolution by integrating the Sentinel-2, Landsat 8, and MODIS satellites. Therefore, combining multi-source remote sensing data provides opportunities for global rice mapping in the future. In addition, the increasing number of Planet satellite images at higher resolutions (3–5 m) could further improve the accuracy of the paddy rice map in the future (https://www.planet.com/, last access: 5 November 2021). To further improve the accuracy of paddy rice map products, more accurate information on cropland and forest masks and crop calendars will need to be developed in the future.

5 Data availability

The datasets of the paddy rice maps for SE and NE Asia from 2017 to 2019 are available on a public repository. A small example of the data can be found at the following DOI: https://doi.org/10.17632/cnc3tkbwcm.1 (Han et al., 2021b). The full product with 5^∘ per grid cell can be downloaded from the following DOI: https://doi.org/10.5281/zenodo.5645344 (Han et al., 2021a). The spatial reference system of the datasets is EPSG:4326. We encourage users to validate this dataset independently. Please note that some small islands are not classified, and thus data for these areas are not available.

6 Conclusions

We constructed a paddy rice map database for SE and NE Asia for 3 years (2017–2019) at a 10 m spatial resolution (NESEA-Rice10) by integrating MODIS and Sentinel-1 data. The paddy rice planting areas in our database were significantly correlated with those from the official statistics. The distribution of paddy rice in the maps was consistent with existing data products. Additionally, our method reduced the effects of mixed pixels and provided more detailed spatial information than MODIS-based paddy rice maps. We demonstrated that multi-sensor data integration has the advantages of improving the spatial resolution of rice maps and reducing mixed pixels. To summarize, we provided more accurate paddy rice maps at continental scales using the improved method for paddy rice mapping.