High-Resolution Mapping of Paddy Rice Extent and Growth Stages across Peninsular Malaysia Using a Fusion of Sentinel-1 and 2 Time Series Data in Google Earth Engine

: Rice is the staple crop for more than half the world’s population, but there is a lack of high-resolution maps outlining rice areas and their growth stages. Most remote sensing studies map the rice extent; however, in tropical regions, rice is grown throughout the year with variable planting dates and cropping frequency. Thus, mapping rice growth stages is more useful than mapping only the extent. This study addressed this challenge by developing a phenology-based method. The hypothesis was that the unsupervised classiﬁcation (k-means clustering) of Sentinel-1 and 2 time-series data could identify rice ﬁelds and growth stages, because (1) the presence of ﬂooding during transplanting can be identiﬁed by Sentinel-1 VH backscatter; and (2) changes in the canopy of rice ﬁelds during growth stages (vegetative, generative, and ripening phases) up to the point of harvesting can be identiﬁed by Normalized Difference Vegetation Index (NDVI) time series. Using the proposed method, this study mapped rice ﬁeld extent and cropping calendars across Peninsular Malaysia (131,598 km 2 ) on the Google Earth Engine (GEE) platform. The Sentinel-1 and 2 monthly time series data from January 2019 to December 2020 were classiﬁed using k-means clustering to identify areas with similar phenological patterns. This approach resulted in 10-meter resolution maps of rice ﬁeld extent, intensity, and cropping calendars. Validation using very high-resolution street view images from Google Earth showed that the predicted map had an overall accuracy of 95.95%, with a kappa coefﬁcient of 0.92. In addition, the predicted crop calendars agreed well with the local government’s granary data. The results show that the proposed phenology-based method is cost-effective and can accurately map rice ﬁelds and growth stages over large areas. The information will be helpful in measuring the achievement of self-sufﬁciency in rice production and estimates of methane emissions from rice cultivation. B.M., R.M.S., N.C.S. and S.G.E.G.; methodology, R. and B.M.; software, F. and R.; validation, F. and R.; formal analysis, F. and R.; investigation, F. and R.; resources, F. and R.; data curation, F. and R.; writing—original draft preparation, F. and R.; writing—review and editing, B.M.; visualization, F. All authors have agreed to the of the


Introduction
Rice is consumed by more than half of the world's population and rice fields account for over 10% of the global cropland area. Globally, rice is the largest plant-based emitter of greenhouse gases [1]. Accurate and up-to-date information on how much rice can be produced is crucial to achieving food, climate, and water security. This is especially crucial in Malaysia, where rice is the main staple food. To ensure national food security, the government of Malaysia has set a goal of a 70% rice self-sufficiency level (SSL) [2]. Measuring the achievement of the SSL rice target requires accurate and up-to-date information on how many rice fields can be cultivated and how much of them can be harvested. Currently, information on rice fields relies on field survey data that are incomplete, time-consuming, and lack spatial details.
However, most remote sensing studies have only mapped the rice field extent. In Southeast Asia, rice is grown throughout the year, and thus, mapping the growth stages of rice fields would be more valuable. Another challenge in the tropics is that the period and frequency of cropping are variable. Each region has its own cropping calendar-for example, East Java, Indonesia, has four cropping patterns, while West Java has three patterns [39]. Thus, this study aims to fill this research gap by applying a phenologicalbased approach to map rice growth stages using a combination of Sentinel-1 and Sentinel-2 at a national scale. The hypothesis is that the unsupervised classification of Sentinel-1 and 2 time-series data can identify rice fields and growth stages, because (1) the presence of flooding during the transplanting can be identified by Sentinel-1 VH backscatter; and (2) changes in the canopy of rice fields during growth stages (vegetative, generative, and ripening phases) up to the point of harvesting can be identified by Normalized Difference Vegetation Index (NDVI) time series.
There are advantages and limitations for each satellite; Sentinel-1, a typical SAR, is not affected by clouds or illumination conditions. However, it has noise (speckles) that makes the map uncertain. On the other hand, Sentinel-2, an optical satellite image, is limited by cloud cover, but it provides a vegetation index directly related to rice growth stages.
Mapping rice fields using remote sensing can be broadly grouped into machine learning techniques and phenology-based approaches. The machine learning approach is based on the calibration of field observations; while accurate, it requires lots of data. The phenology-based approach relies on the knowledge of the rice growth stages, which are related to the time series of reflectance data. In particular, rice fields can be differentiated from other crops due to the presence of standing water during crop transplanting. Based on limited ground-based data, phenology-based classification can detect plant growth stages based on their vegetative cover changes (i.e., [22,40]). Most studies that use the phenologybased method are based on spectral indices and classify such indices using supervised or manual classification. There is still a lack of studies that employ unsupervised classification to differentiate spectral signatures into phenological classes.
The Google Earth Engine platform (GEE, www.earthengine.google.com, accessed on 9 April 2022) is one of the most promising developments for earth science data access and analysis [41], making the processing of multi-temporal satellite data more streamlined. The GEE provides a variety of raster data sets such as Sentinel-1 and -2, Landsat 4, 5, 7, and 8, MODIS, etc. In addition, the new datasets are periodically ingested. Furthermore, the GEE offers an interactive platform for geospatial and machine learning algorithms. Several studies have utilized the GEE platform for rice field mapping [22,29,30,34,35,40].
In a recent study, NESEA-Rice10 [40] produced a rice field extent map covering Southeast Asia (including Malaysia) and North-East Asia. The map was a new achievement because it covered a large area with high resolution (10 m). However, the map still has not been validated in tropical regions. Furthermore, NESEA-Rice10 only provides static information on where rice is grown.
Based on the research gap identified in the literature, this study introduced a phenologybased approach to map rice field extent and the cropping calendar in Peninsular Malaysia, covering over 13 million ha. This was achieved using an integration of Sentinel-1 and 2 time series data using unsupervised classification in the GEE platform. The outcomes of this study will provide an efficient framework for mapping rice fields and their growth stages, which will be valuable in tackling food security issues and mitigating methane gas emissions from rice cultivation.

Study Area
Peninsular Malaysia-also known as West Malaysia-covers an area of 130,598 km 2 , comprising 11 of 13 states of Malaysia (1 • N-6.50 • N and 99.50 • E-104.50 • E; Figure 1). The topography of the area ranges between 0 and 2187 m. The climate is classified as tropical rainforest (Af) according to the Köppen climate classification, with uniform high-humidity hot temperatures and moderate rainfall throughout the year. The average annual rainfall ranges between 1933 mm and 3080 mm [42]. The daily temperature ranges from 23-32 • C and humidity is constantly high (exceeding 68% throughout the year) [43]. Peninsular Malaysia experiences two monsoons, the northeast monsoon (NEM; November to March) and the southwest monsoon (SWM; May to September), and two short transitional periods from April to May and October to November [43,44]. Heavy rainfall events usually occur during the NEM. Maximum rainfall usually occurs during November and December, although substantial rainfall also occurs in the transitional periods (April and October). June and July during the SWM are the driest periods [44,45].

Rice Cultivation Practices
Rice fields in Peninsular Malaysia are irrigated and have two cropping seasons. Crops during the wet season (August to February) are called the main season or are locally known as "tanaman musim utama". Crops planted during the dry season (March to July) are called the offseason or are locally known as "tanaman luar musim" [22,46]. Due to the availability of irrigation, rice fields in Peninsular Malaysia have a diverse cropping schedule throughout the year.

Satelite Imagery Datasets
Sentinel-1 and 2 time-series data were retrieved from January 2019 to December 2020 from the Google Earth Engine (GEE) catalogue. Table 1 presents the characteristics of the Sentinel-1 and 2 data used in this study. Sentinel-1 ground-range-detected (S1_GRD) images with an interferometric wide swath (IW) instrument mode were used in this study. The IW mode images were provided in dual polarization with vertical transmit, vertical receive (VV) and vertical transmit, horizontal receive (VH). Only VH-polarized backscatter data were used in this study because they are more sensitive to rice fields than VV backscatter [22,50].
We only used Sentinel-1A data, as Sentinel-1B data were not available in the study area. The spatial resolution of the images was 10 m. These images were pre-processed with the Sentinel-1 Toolbox, employing thermal noise removal, radiometric calibration, and terrain correction using the SRTM 30 DEM. After pre-processing, each scene was calibrated and ortho-corrected [51].
Sentinel-2 Level-1C (Top-of-Atmosphere Reflectance, TOA) was used instead of Sentinel-2 Level-2A (Surface Reflectance, SR) because artefacts appear in Sentinel-2 Level-2A SR. In addition, Sentinel-2 Level-1C (Top-of-Atmosphere Reflectance, TOA) have been successfully used in numerous studies for surface classification, such as [28,30,32,34,52,53]. Band 4 (Red) and band 8 (NIR) from Sentinel-2 multispectral image (MSI) level 1C data were used to calculate the normalized difference vegetation index (NDVI). The spatial resolution of the images was 10 m. A total of 883 scenes from Sentinel-1A and 3,825 scenes from Sentinel-2A and B were accessed from the GEE.

Generating Regions of Interest (ROI)
First, we generated the regions of interest (ROI) for each rice field location. According to the prior knowledge and information available, which were checked using very high-resolution images and street view photos from Google Earth, 31 rice fields ROIs in Peninsular Malaysia were established. The area of the ROI varied between 100 to 316,000 ha. The distribution of the ROIs can be found at the following link: https://code.earthengine.google.com/8208e211dde42d88108fbe3f5c2c420f (accessed on 9 April 2022).

Monthly Composite Imagery
We then compiled monthly composite data from Sentinel-1 (24 composites) and Sentinel-2 (24 Composites). Sentinel-1 SAR data were not affected by clouds; however, the data had a high noise-to-signal ratio. We used the Median Value Composite (MedVC) of the monthly time series to reduce the noise in the VH polarization data, especially in the area with overlapping data scenes [22,54].
We calculated the NDVI from the radiance on the red (band 4) and near-infrared (NIR; band 8) measured by Sentinel-2 as follows: However, several factors can cause a considerable change to NDVI values. The most important factors are changing illumination and viewing conditions within a single image, and also between images of different days; the presence of cloud cover; and variation in atmospheric constituents-in particular, the variation in water vapour and aerosol concentrations [55]. Maximum Value Composites (MaxVCs) are generally accepted as an effective way to reduce these undesirable factors [55] and have been successfully used in studies using MODIS and SPOT data [56][57][58][59]. Thus, the MaxVCs of NDVI were used to generate the monthly composite imagery from Sentinel-2.
The monthly composite for both Sentinel-1 and 2 also addressed the issue of time alignment as the two satellites acquired images at different times. Using monthly values also allowed alignment with the rice growth stages described in Section 2.2.

Clustering Time Series of VH Polarization and NDVI
We used an unsupervised classification approach to group the time-series data into rice groups [60]. A total of 48 monthly composites of VH polarization and NDVI data were stacked as input into the K-Means model (i.e., the ee.Clusterer.wekaKMeans function in the GEE). The algorithm normalized the numerical inputs and used the Euclidean distance to measure distances between clusters, and iteratively modified the centroids (cluster means) to minimize within-cluster variation (distances) [22]. In total, 3000 random points for each of the ROIs were sampled as training data for training the K-Means model with 20 clusters. The GEE uses a random-point sampling approach to determine the centroids and generate the clusters. Since these points were randomly distributed across the whole region of the targeted ROI, we assumed that they adequately represented the variations of the spectra data. Moreover, the random sampling method does not require more memory computation in the GEE than the stratified random method. These random points and clusters were able to differentiate all land covers in the area (rice field, waterbody, built-up, tree).

Extracting Representative VH Polarization and NDVI Cluster Profiles
To extract representative VH polarization and NDVI cluster units identified by the K-means clusters, a total of 1000 samples were generated randomly for each cluster. Subsequently, spatial median values of the monthly time series for each cluster were used as the representative VH polarization and NDVI time-series profiles.
The time required by the GEE for clustering and cluster profile extraction depended on the area and the speed of internet access. For example, for a ROI of 316,000 ha, it took about 40 s to process the data using a personal computer with 16 GB RAM and a core i7 8th gen.

Labelling and Classification
Rice field and non-rice field areas were differentiated visually based on the representative profiles of VH and NDVI time-series data. Figure 3 shows an example of the temporal evolution of the VH polarization and NDVI cluster profile for four land-cover classes: built-up, tree, water body, and rice field. Rice fields have a unique time-series cluster profile in which the VH polarization and NDVI values fluctuate seasonally. In contrast, signals on other land covers were relatively constant, as shown in Figure 3. The median VH polarization and NDVI values for water bodies were −28.714 dB and 0.040, trees about −15.049 dB and 0.678, and built-up areas about −15.899 dB and 0.185. These values are in agreement with past studies [11,[13][14][15]22,26,61]. The time that was required for interpretive cluster labelling took longer-about half an hour for each ROI-and depended on the number of rice field clusters.

Identifying Rice Cropping Calendar
We determined the rice cropping calendar based on representative profiles of the VH and NDVI time-series data. The rice cropping calendar can be derived from rice cropping intensity and growth stages, including soil tillage and planting, as well as vegetative, reproductive, and maturity stages. For growth stages, first, the soil tillage and planting stage should be identified via VH backscatter and then followed by the next growth stages.

Accuracy Assessment
The accuracy of the produced maps in this study was evaluated using three types of assessment: (1) validation using virtual observations based on Google Earth VHR (very high resolution) and street view images, (2) comparison with existing rice map by NESEA-Rice10 [40], and (3) comparison with known planting and harvesting schedule.
We generated 760 random points across all ROIs across Peninsular Malaysia for the first assessment. After this, we manually evaluated each location based on VHR images from the year 2020 and Google Earth street view images, where 348 points were identified as rice fields and 392 points as non-rice. The accuracy of the rice maps was then evaluated in these virtual field data using the confusion matrix and Kappa value.  In the second assessment, this study's overall accuracy and kappa coefficient were compared with those from the NESEA Rice 10 map. Finally, the third assessment compared the cropping calendar generated from this study with reference data from the granary area in Selangor state (IADA BLS) and Kedah state (MADA).

Spectra Characteristics of Paddy Fields
Based on rice phenology, the response of the Sentinel-1 and 2 spectra to three key phenological periods (Figure 3) can be described as follows: 1 The tillage and planting (T) period, where rice fields are flooded. According to the VH polarization data in Figure 3, the rice fields had a local minimum VH polarization value about <−24.5 (i.e., close to the value of a water body). This value depends on water levels or soil moisture during soil tillage. Higher water levels lead to more negative VH values [22]. On the other hand, the NDVI value during this period may not be at the required minimum due to cloud cover. Thus, Sentinel-1 data (VH polarization value) are more sensitive to detect the transplanting period than Sentinel-2 data. 2 The vegetative (V) and reproductive (R) period, where rice grows rapidly with canopy closure. Consequently, increasing chlorophyll signals are accompanied by a significant decrease in soil signals [29]. NDVI values rapidly increase during this period, peaking at the end of the generative period, when the crop enters the heading or maturity stage (NDVI value was about 0.7). On the other hand, VH polarization also increases in this period and continues to rise until the end of phenology. 3 The maturity (M) period. Characterized by a rapid decrease in chlorophyll content as the carotenoid content increases [29]. This distinctive drop in NDVI values indicates the maturity of the rice plant. This period occurs in the last month of the rice growth cycle.
Sentinel-1 data (VH polarization value) were more sensitive to detecting the transplanting period; on the other hand, Sentinel-2 data (NDVI value) were more detailed in showing the stages of growth and maturity. The difference between the local minimum and the peak was about four months, representing one growing season. The VH polarization and NDVI values decreased during the fallow period. Figure 3 shows an example of a rice field cluster in Selangor that had a double cropping season in a year, where the first (October-January) and second (April-July) planting seasons were referred to as the main and off planting season, respectively. The time series profiles of VH polarization and NDVI can be designated with the four periods: T: tillage and planting (30 days); V: vegetative (30 days); R: reproductive (30 days); and M: maturity (30 days). Thus, one season of rice cropping was about 120 days or 4 months.

Accuracy of Rice Maps
Based on the unsupervised classification, rice areas throughout 31 ROIs across Peninsular Malaysia were identified. Maps of the rice extent in Peninsular Malaysia from 2019-2020 are presented in Figure 5. The complete map is available at https://rudiyanto. users.earthengine.app/view/malaysiarice, accessed on 9 April 2022. The maps produced in this study successfully distinguished rice fields from other land covers.   We validated our map using random sampling in the ROI of rice fields to evaluate the accuracy at a pixel level. Our map had an excellent overall accuracy of 95.95% and a kappa coefficient of 0.92, as shown in Table 2. We also evaluated the accuracy of the NESEA-Rice10 product on the same random samples. NESEA-Rice10 had a lower overall accuracy of 90.27% and a kappa coefficient of 0.80 (Table 3). Omission errors occurred when rice fields were classified as non-rice fields. Commission errors occurred when non-rice fields were classified as rice fields [22]. This study's omission and commission errors of rice fields were 1.1% and 7.0%, respectively, while NESEA-Rice10 errors were much larger-19.9% and 1.1%, respectively. We note that outside the ROI, the NESEA-Rice10 had many speckles where non-rice areas were classified as rice.

Cropping Calendar
A total of forty different cropping calendars of rice fields were identified across Peninsular Malaysia based on K-means spectra profiles of NDVI and VH polarization. The cropping phenological parameters of growth stages for each rice field group from January 2019 to December 2020 are listed in Table S1. Figure 6 shows the representative maps of cropping patterns with information on rice groups in each state corresponding to the cropping pattern in Table S1. For example, Selangor state had four rice groups, Kedah state had two rice groups, Perlis state had one rice group, Penang had two rice groups, and Besut district, Terengganu state had one rice group.   Table S1. Table 4 shows an example comparison of the cropping calendar from the published granary area schedule and this study. The cropping calendar derived from this study agreed well with the granary area schedule. Peninsular Malaysia has two cropping seasons, the main season and off season, with no cash crops in-between seasons. As an example, Selangor and Kedah states had two cropping seasons in one year. In Selangor, there were four cultivated areas with different cropping schedules. The main season in 2020 in Selangor started in July, August, September, and October for Group 2, 1, 4 and 3, respectively. Concurrently, the off season started in January, February, March and April, respectively. In Kedah, two cultivated areas had different cropping schedules; the main season of both rice groups started in October, while the off season started in April and May for Groups 2 and 1 ( Table 4). The Kedah and Selangor state cropping calendar correlated with the rice variety planted in the area, MR 220 CL2 (97-113 growing days).
From the rice field extent and cropping calendar data in Table S1, an aggregated number of planted and harvested rice field areas in Peninsular Malaysia from 2019-2020 were generated and are shown in Figure 7. In 2019, the largest rice area was harvested in February at 143,980 ha, while in 2020 the largest area was harvested in August at 155,481 ha. There was an anomaly in April 2020 caused by the Kedah-2 cluster, which had 5 months of growth. The planting period was in April and the cluster was harvested in August over 54,782 ha. This shift was caused by rice fields that were planted at the end of April 2020 and the beginning of May 2020 and were harvested at the same time in August 2020. In 2019, Kedah-2 had 4 months of growth that started in May and finished in August. As a result, the harvested rice area in 2019 was 458,064 ha, while in 2020, it increased to 460,388 ha. Table 4. Comparison of cropping calendars from granary area schedules and this study. Table 4 corresponds to Figure 6a,b. Table 4 corresponds to Table S1. T = Tillage, V = Vegetative-1, V2 = Vegetative-2, R = reproductive, M = Maturity, F (blank) = Fallow. Orange colour is the main season and green colour is the off season.   Figure 7 corresponds to Table S1.

Targeted ROI Processing Approach
Mapping rice fields at a fine resolution is challenging when applied to large areas such as Peninsular Malaysia. A targeted ROI processing approach was applied to reduce computational costs. This approach also allows the running of the algorithm in parallel. Rice field locations can be obtained from agricultural statistics such as those of the local government [46] or international organizations-for example, rice maps in Asia from the International Rice Research Institute (IRRI) [40]. As the rice-growing areas in Malaysia are known, the ROI can be generated for every known area. This approach significantly reduced the processing time as the classification process focused on the known rice areas to pinpoint and refine rice fields.
The targeted ROI processing approach also reduced commission errors and "salt and pepper" effects. Conventionally, non-rice areas are masked using heuristic rules to reduce the salt and pepper effect. For example, masking areas with a slope greater than 5 • to remove steep terrain areas that are unsuitable for rice [40,62]; masking forest lands using Global PALSAR Forest Map 2017 [63,64]; excluding wetlands because paddy rice fields and wetlands have similar characteristics in flooding signals [65,66]; or masking sparse and natural vegetation using annual maximum vegetation index values that are less than 0.3 [35]. Although non-rice areas were reduced, other land covers may still be included in the data processing. As a result, the resulting maps still had a large amount of noise, as demonstrated by a comparison between the rice map of this study and that of NESEA-Rice10 [40] in Selangor State, Malaysia (Figure 8d). The map produced in this study shows a clearer separability between rice and non-rice areas such as roads, irrigation channels, built-up areas, and other vegetation covers (Figure 8a,b). The salt and pepper effects were significantly reduced (Figure 8c). We can see that in NESEA-Rice10, some waterbody and other vegetation covers were misclassified as rice fields (Figure 8d), which led to an overestimation of rice field extent. The targeted ROI processing approach also reduced the image processing time. Since K-means clustering requires many sampling points, the random sampling points in the targeted ROI approach could pick up small rice fields. This approach had an advantage in the classification of fragmented rice fields. However, accurate rice field regions need to be known beforehand in order to implement ROI-targeted processing.

Combination of Sentinel-1 and 2
This study used a combination of Sentinel-1 and 2 time series to capture the rice extent and the dynamics of rice growth. Tropical regions tend to have many clouds that interfere with optical satellite images. The Maximum Value Composite (MaxVC) of NDVI is generally used to minimize the effect of cloud cover. In rice mapping studies, the MaxVC of the NDVI has successfully minimized the effect of cloud cover in MODIS data [56][57][58]67], SPOT vegetation data [59], and Landsat-8 and Sentinel-2 data [28,34,68]. As the NDVI of water is smaller than the NDVI of clouds, cloud cover could be mistaken for water. Thus, the MaxVC of the NDVI cannot effectively detect rice fields during the transplanting period, as shown in Figure 3b. A limitation of the MaxVC of the NDVI is when there is persistent cloud cover over the study area during all revisiting times (a 5-day period) within a month. To overcome this problem, spatial interpolation [69] of NDVI values in the regions affected by clouds could be performed using machine learning models using SAR imagery and other environmental covariates [70].
On the other hand, the median of the VH polarization can detect rice fields during the transplanting period, as shown by the minimum peak in Figure 3a. However, the MaxVC of the NDVI showed better spectra discrimination than VH polarization during other growth stages. Thus, a combination of Sentinel-1 and 2 improved the overall results. The use of monthly greenest images or maximum NDVI value composites of Sentinel-2 data circumvents the need for masking clouds and cloud shadows, as is commonly carried out in rice mapping studies [64,71].

Phenology-Based and Unsupervised Classification Methods
Machine learning and deep learning approaches can be applied as classifiers when mapping rice extent and cropping calendars [10,[72][73][74][75][76]. Since they are supervised classification methods, they require a large number of observations for model training. For example, a recent study [52] used 21,696 observations to train a convolutional neural network (CNN) model to map rice extent and growth stages in West Java, Indonesia. On the other hand, the phenology-based method used in this study does not require many observations. Our method based on k-means clustering can quickly identify rice fields from expert interpretation due to the unique spectra pattern of VH polarization and the NDVI of rice fields from Sentinel-1 and 2, respectively, resulting in high-up to 96%-accuracy rice maps. Due to its ability to generalize clusters from the data pattern, the k-means method can learn from the data and group time series into unique spectral responses related to rice phenology or other land uses. The method was also able to identify reflectance signatures due to soil tillage and planting, vegetative, reproductive, and ripening stages, as shown in Section 3.1. Nevertheless, to identify the reflectance signature, prior knowledge is required. To overcome this limitation, K-means and supervised methods such as random forest [75] can be combined to automatically produce rice extent and growth stage maps.
An important feature of this study was the ability to map cropping intensity and cropping calendar data. RiceAtlas [77] and RICA [78] provide rice calendar and production data from around the world and rice crop calendars for Asia, respectively. However, those studies only present information on cropping calendars at the country level. In reality, the cropping calendar varies significantly within countries and could change due to local conditions. Using a phenology-based method, we have provided not only maps of rice extent but also cropping intensity and calendar at the pixel level. We have demonstrated that the unsupervised classification of Sentinel-1 and 2 time-series data and their implementation in the GEE is cost-effective in mapping rice extent, intensity, and calendars. In addition, the approach is scalable, and GEE can process data of vast areas [40,79].

Mapping Rice Extent in Tropical Regions
Our methods were applied in rice fields in tropical regions that are planted throughout the year with different planting dates and frequencies. In South-East Asia, different places can have different rice cropping calendars due to variable irrigation schemes. The novelty of our method is that it can handle the complexity of rice fields in tropical regions. Most previous studies have used Sentinel-1 and 2 for mapping rice in temperate regions such as Japan [30], China [32][33][34][35][36][37][38], and India [31], which are only planted in summer. Our methods can map rice fields rapidly-in less than an hour for 100,000 ha-as rice fields can be easily identified from their unique phenology spectra without extensive surveys. Thus, our method has a high potential and capability for application in other regions.

Sources of Uncertainty
Although the present results clearly show the high accuracy of the rice-mapping procedure, it is appropriate to recognize several potential limitations. Errors usually occurred at the perimeter between rice and non-rice fields (Figure 9), leading to omission errors in the rice field category. The mixed pixel effect in the boundaries along the field edge contributed to the discrepancy. This effect has been reported in rice mapping studies and is still a source of spatial uncertainty [3,15,20,22,50,73]. For the commission errors in the rice field category, the contribution was mainly from non-rice fields (such as small roads, levees, or bunds) classified as rice fields. Bunds are similar to rice fields, as their water levels and soil moisture also fluctuate seasonally.

Conclusions
Our results support the hypothesis that the integration of Sentinel-1 and 2 time-series data using a phenology-based method in the GEE cloud computing platform can accurately produce maps of rice extent, cropping intensity, and cropping calendars. Furthermore, the proposed method produced rice extent and growth-stage maps with a high spatial resolution (10 m) over a large region.
The results show that the proposed method had a high overall map accuracy of 95.95%, with a kappa coefficient of 0.92. The method distinguished rice growing stages and established that most of the rice fields in Peninsular Malaysia had two cropping seasons a year.
The proposed methodology had several advantages compared to past studies: • It produces not only high-resolution maps of rice extent but also of intensity and cropping calendars in tropical regions due to the use of both Sentinel-1 and 2 time series data; • It produces high-accuracy maps rapidly without extensive field survey data due to the use of an unsupervised K-Means method in the GEE; • It reduces salt and pepper effects due to the use of targeted regions of interest; • It eliminates the need for masking clouds and cloud shadows due to the use of monthly greenest images or maximum NDVI value composites of Sentinel-2 data; • It is able to identify the signatures of rice phenology stages, including soil tillage and planting, vegetative, reproductive, and ripening stages due to the fusion of Sentinel-1 and 2 time series data.
The results presented in this paper could be used as a complementary source of information to field survey data. This information is necessary to achieve the United Nations 2030 Agenda for the Sustainable Development Goals (SDGs) [80]-particularly Goal 2, Zero Hunger. As paddy rice consumes a large amount of water, this data can be used to address water security management (SDG Goal 6). Rice fields are also an important source of methane, and the spatial information can be used for greenhouse gas accounting (SDG Goal 6). Finally, the data can help form government policies to address national food security to reduce poverty (Goal 1).