Phenology-based delineation of irrigated and rain-fed paddy fields with Sentinel-2 imagery in Google Earth Engine

ABSTRACT Paddy rice agriculture is practiced in both rain-fed and irrigated ecosystems in the Philippines. However, small farms are prevalent in the region, and current satellite-based mapping techniques do not distinguish between the two ecosystems at farm scales. This study developed an approach to rapidly map irrigated and rain-fed paddy rice in Iloilo, Philippines at 10 m resolutions using Google Earth Engine. This approach used an ensemble of classifiers based on time-series vegetation indices to produce dry and wet seasonal maps for the entire province. Results showed a predominance of rain-fed rice areas in both seasons, with irrigated rice making up only one-fourth of the total rice area. The overall accuracy was achieved at 68% for the dry season and 75% for the wet season based on ground-acquired points and very high-resolution imagery. The two types of paddies were classified at accuracies up to 87%. Furthermore, the land cover maps showed a strong agreement with the municipal statistics. The resultant maps complement current official statistics and demonstrate the prowess of phenology-based mapping to create paddy inventories in a timely manner to inform food security and agricultural policies.


Introduction
The world population is projected to increase to 9.77 billion by 2050 (UN-DESA 2017). A larger population will impose a higher demand for food, particularly staples such as rice (Oryza sativa L.). Rice is the world's most important food crop and is primarily grown for human consumption. It is cultivated chiefly in the lowlands of developing nations in Asia (GRiSP 2013). Rice yield is affected by many factors, including seasonal and annual climate events such as flooding, typhoons, frost, dry spells, and drought, along with soil erosion and pestilence (IPCC 2013(IPCC , 2018Wassmann et al. 2009). Besides, rice yield is also affected by the manner of cultivations, such as rain-fed or irrigated. Thus, the accurate rice yield estimation requires the mapping of rice crops into rain-fed and irrigated classes. Monitoring the seasonal changes to areas devoted to different rice ecosystems is integral to ensure food security and sustainable development in years to come.
Crop yield statistics are routinely collected at various government levels but have long-standing limitations. In the developing world, subjective methods such as eye estimation, expert assessment, and farmer interviews are used (FAO 2017), but these traditional methods are imprecise. Ground-based methods utilizing crop cuts or Global Positioning System-enabled devices can improve the quality of estimates, but costs and resources are significant limitations when data is routinely needed (FAO 2016;Kabir et al. 2016).
Space-borne satellite remote sensing addresses these limitations, offering comprehensive area coverage with frequent repeat observations. Imagery from these sensors provides a plethora of data detecting biophysical characteristics of rice such as biomass, water content, and photosynthetic activity (Kuenzer and Knauer 2013). Several rice-focused studies have used optical sensors such as SPOT-VGT, MODIS, AVHRR, and Landsat (Manjunath et al. 2015;Zhang et al. 2015;Guan et al. 2018;Thenkabail et al. 2007;Bachelet 1994). These sensors feature large swath widths, frequent observations, and long-term data records, which are ideal for crop mapping. However, the spatial resolution of these sensors range from 30 m to 1 km. This range is too coarse to effectively detect spatially complex rice environments, such as smaller farms that predominate rain-fed ecosystems (Koirala et al. 2014;Palis 2020;Adamopoulos and Restuccia 2014). Furthermore, mixed pixels cause mapping accuracy to decline in coarse-resolution sensors as they contain spectral signals of multiple covers other than paddy fields (Okamoto and Fukuhara 1996). One solution is to use finer resolution images, such as WorldView and RapidEye (Wan and Chang 2019;Kim and Yeom 2014). However, these are costly to acquire and have longer gaps between repeated observations. Sensors with high temporal resolution allow accurate mapping of distinctive temporal profiles of rice and its flooding patterns (Liang et al. 2019;Zhang et al. 2015;Gumma et al. 2014). Recently, constellations of satellites sought to fill this gap, notably the Copernicus Program's twin Sentinel satellites. These have 10 m resolutions and 5-day repeat intervals. Several studies have demonstrated the suitability of these sensors for rice mapping (Ferrant et al. 2017;Jin et al. 2017;Zhang et al. 2018;Son et al. 2020).
Mapping methods using dense time-series imagery are standard for analyzing rice distribution and yield (Shew and Ghosh 2019;Mosleh, Hassan, and Chowdhury 2015). Two problems persist with singledate optical imagery of rice areas: (1) frequent cloud cover and (2) highly dynamic plant cover throughout the growing season. Multi-date imagery solves this by tracking critical growth stages, while statistical techniques can synthesize cloud-free datasets (Zhang et al. 2015;Jin et al. 2016;Teluguntla et al. 2015). This approach exploits the unique phenology of rice. Rice is essentially a kind of vegetation, so it is no different from other vegetation if detected from a vegetation index. However, rice is unique because it has a short growth period, with distinctive stages of growth at which its spectral signature follows a definite pattern different from that of other vegetation. Most studies utilize the Normalized Difference Vegetation Index (NDVI), but it saturates at denser canopies with moderate-to-high biomass (Gitelson 2004). Saturation affects the ability to detect temporal changes and monitor vegetation dynamics, leading to confusion of covers with similar growth patterns.
Researchers developed alternative indices to be more sensitive at higher biomass stages. The Green Chlorophyll Vegetation Index (GCVI) -was initially designed to separate soybeans from maize and was linearly correlated with chlorophyll concentration (Gitelson et al. 2005;Baret and Guyot 1991;Gitelson, Gritz, and Merzlyak 2003). GCVI was used to map irrigated croplands across the continental US, resulting in low mapping errors and bias (Ozdogan and Gutman 2008). Despite its simplicity and advantages, current literature rarely used it to map rice areas despite its superior ability to detect canopy chlorophyll and the prevalence of rain-fed/irrigated ecosystems in rice agriculture (Guan et al. 2018;Zhang et al. 2018). Another approach to rice mapping combines a vegetation index (VI) with a water index such as the Land Surface Water Index (LSWI) (Gao 1996;Xiao et al. 2002). A modified form of NDVI, this method effectively detects the flooding-transplanting signal of rice paddies. LSWI and NDVI have been jointly used for phenology-based classification of rice in several regions of Asia, albeit with varying success (Jin et al. 2016;Xiao et al. 2005Xiao et al. , 2006. Over large areas, cloud-based platforms such as Google Earth Engine (GEE) have gained increased use as an effective platform for analyzing time-series data due to its vast catalog of analysis-ready geospatial datasets combined with high-performance computing capabilities (Gorelick et al. 2017). Several studies demonstrated GEE's ability to map broad land covers, such as crop/non-crop or rice/non-rice areas, but none attempted to detect irrigated and rainfed rice in tropical regions using high-resolution imagery (Singha et al. 2019;Clauss et al. 2018;Lee et al. 1994;Guan et al. 2018;Zhang et al. 2018;Teluguntla et al. 2018;Oliphant et al. 2019). Separating these two ecosystems is essential as they have different crop and water management strategies. In particular, rainfed farms are particularly vulnerable to climate change, tend to have lower yields and smaller plot sizes (Koirala et al. 2014;Palis 2020;Adamopoulos and Restuccia 2014).
Successful mapping of rice ecosystems depends on the use of a suitable classification approach. The predominant method is using single classifiers such as a supervised or unsupervised classifier . However, traditional classifiers designed for spectral information may not be suitable for the phenology-based temporal approach to a complex environment such as tropical paddies. A solution to this is to combine the strengths of several classifiers at several stages for the mapping (Du et al. 2012). Additionally, knowledge-based approaches using information external to imagery such as irrigation canal networks or crop calendars can improve the mapped results (Ragettli, Herberz, and Siegfried 2018). Nevertheless, their effectiveness in the mapping remains unexplored. There is a gap in the current literature in providing a timely and accurate differentiation of the spatial extent of rain-fed and irrigated rice agriculture at a comprehensive areal coverage with a resolution that can detect smallholder farms. Therefore, the overarching objective of this study was to investigate how effectively the use of spectral and temporal knowledge gathered from satellite imagery and external knowledge can improve the mapping of irrigated and rain-fed rice in Iloilo Province, the Philippines. Specifically, our aims were to (1) investigate the vegetation index temporal profiles for different land covers in the study area, (2) map the two main rice ecosystems and quantify their areas and changes between seasons, and (3) assess the accuracy of the maps with field data and official statistics. Although this study was conducted in a single region and year, this time-and cost-effective methodology is envisioned to be adapted to more expansive areas and extended to multiple years.

Study area
The study area is Iloilo, the largest of four provinces located on the island of Panay in the Western Visayas group of islands of the Philippines (Figure 1). It extends from 10.47 • N to 11.65 • N and 122.01 • E to 123.38 • E, covering a total land area of 5,079 km 2 (Iloilo PPDO. 2019). Its topography consists of flat plains, mountainous ranges, and rolling hills. Onethird of the province is level plains are full of fertile agricultural land with extensive river systems. To the west lies the Central Panay Mountain Range. The highest point is Mount Baloy at 1,908 m. Loam is the predominant soil type in the province and is considered suitable for most crop types (Parida et al. 2008;Iloilo PPDO. 2019). Iloilo has a tropical monsoonal climate with two distinct seasons. The dry season occurs December-April, and the wet season lasts for the remaining months. Average monthly temperature ranges from 25.8 to 28.5 • C. Monthly precipitation ranges from a low of 27 mm in February to 346 mm in August (Iloilo PPDO. 2019).
The agricultural landscape of Iloilo is predominantly paddy rice, but other crops are present. Corn is found mainly on the hillsides or as a second crop in some areas, and sugarcane, primarily in the central municipalities. Perennial crops such as banana, coconut, and coffee are also prevalent in mountainous regions. Major paddy areas are located in the east and central sections of the province. Rain-fed and irrigated rice are the main rice ecosystems in Iloilo, with some upland rice. Rainfed rice is dependent on seasonal rains and is typically planted once a year. Irrigation infrastructure is less prevalent in Iloilo, with most concentrated near river systems, such as the Jalaur River, but there are also small-scale or private irrigation systems. Areas with irrigation can plant rice twice a year. Upland rice paddies with smaller plot sizes dominate the southern and western regions near the main mountain ranges. The main cropping season starts in June and lasts until the September harvests. The second cropping season typically lasts from October to January. In 2017, the province produced 937,268.72 MT of milled rice, 39% of which comprised wet season rain-fed rice, totaling 363,899 MT (Philippine Statistics Authority 2018). Evergreen forests are primarily found in the mountainous regions, while coastal areas to the west are lined with mangrove forests and fishponds. Iloilo City is the primary urban center in the central-west portion of the province.

Data used
The Sentinel-2 Level 1C data used in this study are the top-of-atmosphere (TOA) reflectance product, available from the GEE Data Catalog as "COPERNICUS/ S2". From digital numbers, the image is radiometrically corrected for radiances at TOA values. Corrections applied were dark signals, pixel response non-uniformity, crosstalk, interpolation of defective pixels, restoration of high spatial resolution bands, and spatial filtering for 60 m bands. The geometric correction was then applied using a Digital Elevation Model (DEM) to project imagery in cartographic coordinates. The imagery was then associated with pre-defined 100 km 2 x 100 km 2 tiles and resampled to the native resolution of each band (orthorectification). Radiances were converted to TOA reflectance values at this point.
To coincide with the main growing seasons of rice of 2019 in the region and the timing of official harvest statistics surveys, we used data from October 1, 2018 to October 1, 2019. Four Sentinel-2A/B scenes covered the whole province. This number rose to 144 over the study period ( Figure 2). All the images were cloud-masked using the in-built cloud mask band of Sentinel-2 Level 1C products, Figure 1. Location of the study area, the province of Iloilo. Names of municipalities of the province are given within the administrative boundaries, delineated by red borders. Elevation data shown is from ALOS World 3D (AW3D30) .
QA60. This band contains information on the presence or absence of opaque clouds and cirrus clouds for each pixel. Clear pixels were selected in this instance. The shapefile from the Philippine Statistics Agency (PSA) provided the geographical bounds of Iloilo Province.
For ancillary data, we used the ALOS World 3D (AW3D30). We present a digital surface model (DSM) dataset with a horizontal resolution of 30 m (~1 arcsecond). Spatial data from irrigation networks were acquired from the National Irrigation Administration -Region VI to assist in labeling clusters from classification and validating the results.

Reference data and ground observations
A reference sample site database was created through a combination of fieldwork and high-resolution imagery within GEE. Fieldwork was undertaken throughout the study area from August to November 2019 ( Figure 3). A team of local officials and agricultural technicians visited 43 rice field plots in 18 municipalities across the province and interviewed the farmerowners of each farm. Plots were selected by agricultural officers who had expert knowledge of each farm's representativeness across their municipality. Data collected during the campaigns included: (1) geographic  . Location of reference sites (n = 463) used in this study overlain on a false-color median composite Sentinel-2 image using Bands 2, 8, and 11. Reference sites consisted of (1) site visits with interviews to 43 farms and (2) only geographic coordinates with additional validation using a custom GEE app for 420 sites.
coordinates of the plots using a Garmin GPSMAP 62s handheld unit, (2) rice crop type separated by rain-fed or irrigated rice, (3) crop planting and rotation schedules, (4) farm area measured using GPS and reference fine-resolution imagery, and (5) geotagged photographs of the farm's setting and condition.
An additional set of 420 points was collected across all municipalities in Iloilo using only GPS coordinates, together with information on the land cover type, presence of irrigation infrastructure, and geotagged photographs. The quality of these observations was assured with the assistance of a customized GEE application and Google Earth Pro's Street View imagery. The custom application ( Figure 4) graphed the temporal profile of the candidate location. Additional validation, as needed, was achieved through Street View. Points were excluded if they were misidentified, occupied a small area, or recently changed the land cover. Phenological profiles derived from these samples provided the ideal time series for each class, namely rain-fed rice, irrigated rice, non-rice crops, built-up /bare soil, forest/trees, and water bodies.
Statistics of rice cultivation area for 2019 were acquired from the Provincial Agricultural Office in Iloilo for all municipalities of the province (n = 44). Additional sub-national statistics were obtained from the Palay Production Survey of the Philippine Statistics Authority.

Image processing and classification
Rain-fed and irrigated rice paddies were mapped for the dry and wet cropping seasons of 2019 using a combination of classification methods at 10 m spatial resolution ( Figure 5). Satellite images from Sentinel-2A/B were first merged into 10-day maximum value composites, from which NDVI and LSWI were then computed using Bands 3, 8, 8A, and 11 for the whole year (Table 1). Several temporal and statistical subsets were created. For each seasonal map, unsupervised k-means clustering was applied. After all clusters were labeled into water, forest, and built-up areas, a random forest classifier was used to separate rice from the non-rice crop. Ancillary information, including elevation and irrigation infrastructure, was used to delineate irrigated from rain-fed rice. Finally, the produced rice map was assessed and validated for its accuracy through a constructed confusion matrix and comparison with the validation dataset and the municipal statistics.

Spectral indices and temporal compositing
Two spectral indices were derived from the multispectral bands of Sentinel-2 that take advantage of different optical radiative properties of the land surface. With Sentinel-2, GCVI uses the near-infrared (B8) and green bands (B3) and is calculated as: The land surface water index (LSWI) uses the rededge (B 8A ) and shortwave infrared bands (B 11 ), and its equation is given as: Several temporal composites were constructed from the original data stack, known in GEE as an ImageCollection, to track changes in the derived VIs. Decadal (10-day) maximum value composites (MVCs) were assembled to observe the phenological profile with the least amount of cloud contamination for each reference sample and to account for the variation in planting date. With the assistance of the reference data from both in situ observations and a customized GEE app, the average temporal profile was extracted for each spectral index for non-rice, rain-fed rice, and irrigated rice cover types present inside the study area. Their statistics were then computed using GEE's ee. Reduce methods, including the mean, median, maximum, minimum, and standard deviation. These were primarily for unsupervised classification of separate non-rice types.

Clustering analysis
The WEKA k-Means clustering method within GEE (ee.Clusterer.wekaKMeans) was used to partition the data into subgroups so that the sum of withincluster variations is made as small as possible (James et al. 2013;Arthur and Vassilvitskii 2007). The target number of clusters was set to 30, while other parameters were kept at default values. Input features were statistical composites from temporal data of the two indices. We created two subsets for the dry and wet seasons from the original annual stack that covered the entire study dates. The dry season rice crop covered 6 months: 10-01-2018 to 04-01-2019. The wet season stack ranged from 04-01-2019 to 10-01-2019. Afterward, we generated mean, median, and standard deviation composites of each stack for both spectral indices. Then, we used the percentile reducers within GEE to produce quantile composites at 5%, 25%, 50%, 75%, and 95% for GCVI and NDVI. The resulting multidimensional image array for each stack resulted in a 16-band statistical composite image. For instance, the wet season stack had the following bands: GCVI_mean, GCVI_median,GCVI_stdDev,GCVI_p5,GCVI_p25,GCVI_p50,GCVI_p75,GCVI_p95,LSWI_mean,LSWI_median,LSWI_stdDev,LSWI_p5,LSWI_p25 LSWI_p50,LSWI_p75,LSWI_p95. The final input feature for clustering was a combined wet and dry season stack that had 32 bands.  The annual phenological profiles of GCVI and LSWI were extracted for the resulting clusters. A decision tree was constructed using thresholds determined from one of these indices, primarily GCVI, or its combination with LSWI from the ideal profiles from field surveys. The rulesets could quickly group similar clusters.
However, several clusters did not match the ideal temporal profiles. These were assessed against field validation points collected using the custom GEE application. The cluster map was overlaid on highresolution Google Earth satellite imagery, annual statistical composites, and true-color Sentinel-2 imagery. These were then manually assigned to the corresponding class.
After matching the temporal profile of each cluster with the reference profiles or through manual labeling, the classes were labeled as water, built-up/bare, forest/ trees, and cropland. The output was spatially filtered using the in-built morphological operators within GEE over a 3 × 3 weighted kernel to refine the classification results and reduce the salt-and-pepper effect.

Random forest classification of paddy fields
A machine learning method, random forest, was applied to the cropland cluster produced from the previous step. The algorithm is available within GEE as ee. Classifier.smileRandomForest. The reference dataset was split into 70% training (n = 336) and 30% validation samples (n = 143). Parameters for the classifier wawere00 trees and a minimum leaf population of 10, while others were set to default values. The resulting two classes were non-rice cropland and rice cropland.

Determination of irrigated and rain-fed areas
The rice cropland class was then further classified using a knowledge-based approach. The National Irrigation Administration provided irrigation shapefiles for the main and lateral canals. Elevation was taken from AW3D30. Slope gradients were computed from this dataset using the in-built ee.Terrain.slope algorithm in GEE.
A decision tree was constructed where irrigated rice fields were within a 1 km radius of the irrigation canals, at or below 50 m above sea level, and at a slope ≤3º. These thresholds were based on previous research and also on inspection of the maximum elevation of the canal network dataset (Brouwer et al. 1988;Xiao et al. 2005Xiao et al. , 2006Ghazaryan et al. 2018).

Accuracy assessment
Using the validation subset, confusion matrices and accuracy metrics were computed for the wet and dry season maps of rice/non-rice cropland and rain-fed /irrigated rice paddies. The identity of the selected pixels was then compared to that in the reference dataset, and the results were shown in a confusion matrix, from which four accuracy metrics were calculated, overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA), and kappa coefficient. Additional accuracy indication was obtained by comparing the mapped area of rice areas with the provincial government statistics of the same year. We determine percent difference from official statistics using Equation (3), while the percent area change between dry season and wet season area was computed using Equation (4): (4) where CA is the area of paddies classified from the Sentinel data and OS is the official statistical area of rice harvested from the Philippine Statistics Authority (PSA). WS area is the harvested, wet season area, and DS area is the harvested, dry season area obtained from the images.

Phenological profiles of different classes
GCVI temporal patterns differed from each other among the mapped land covers (Figure 6). Built-up and water bodies had low GCVI values (<0.5) and low monthly variability. For non-crop vegetation areas, the highest GCVI values were from forested areas or trees. For crops and grasslands, seasonal variations were observable. Sugarcane and corn, which have a high biomass, can reach GCVI > 3.0 at peak biomass and consequently decrease 1.5-2.0 after harvest. Grassland has a monthly variation of GCVI values is less acute during the dry season.
Rain-fed and irrigated rice also displays distinct monthly variations. The profiles both showed a double-cropping pattern, but low and high values were different due to various cultivars planted and the timing of planting. For the dry season, the GCVI peaks in December then declines sharply. According to interviews with farmers in the field, these changes coincide with the approximate timing of heading and harvest. Irrigated rice is planted later in December and harvested in the following February/March. For the wet season, rain-fed rice is grown in May/June and harvested in August/September, while irrigated rice starts in July/August. In Iloilo, farmers typically plant 90-day and 100-day varieties, and this is consistent with the length of the growing period observed from the time-series profiles.

Mapped rice paddy
For the dry and wet crop seasons spanning October 2018 to October 2019, our ensemble classification approach identified six land-use types: forest/ trees, urban/bare, water, non-rice crops, rain-fed rice, and irrigated rice (Figure 7). Table 2 summarizes area estimates for rice and non-rice, including changes from dry cropping season to the wet season. Dry season rice had an area of 122,240 ha for rice and 59,051 ha for non-rice crops. Rice cover was composed of 79.07% rain-fed and 20.93% irrigated. The wet season rice area was 103,571 ha and the non-rice crop area of 86,482 ha, of which 75.72% was rainfed, while 24.28% was irrigated. Rice fields account for 25.92% of the province's total land area in the dry season and 21.97% in the wet season. There is a 15.26% decrease in rice area from the dry season to the wet season, while non-rice crops increased by 46.45% during the same period.
Large, contiguous rice areas are located in the central-eastern portion of Iloilo, around the periphery of the main urban area of Iloilo City (Figure 7). Irrigation facilities of the Philippine National Irrigation Administration (NIA) were also located primarily in this region. The north-eastern areas also have large rice areas switching between rice to non-rice crops in the wet season. Rain-fed rice plots are highly fragmented and located closer to areas with higher elevations and steeper slopes (Figure 8). The province's centralwestern and southwestern areas were predominantly forested areas with much smaller plots of croplands close to the coast or near rivers. Figure 9 shows the distribution of rice area per municipality. Municipalities with the most rice areas are located in the central-west to north-western portions of the province, where the majority of the province's plains and the flat regions are. South-eastern and eastern municipalities and a few northernmost municipalities have the least amount of rice areas. These are predominantly forested areas with steeper   slopes, and more non-rice crops are planted here, such as sugarcane and corn. The rice area showed a general decreasing trend from dry to wet season in 2019. Figure 10 shows areas across the province where changes to crop areas were observed. The most considerable changes were the conversion of non-rice crop classes to rain-fed rice during the wet season. To a smaller extent, non-rice also changed to irrigated in the central-eastern area. There was a very low incidence of switching from rice to non-rice crops in the same period. Table 3 shows the confusion matrix of the mapped rain-fed rice, irrigated rice, and non-rice crop classes in the two seasons. The dry season crop had an overall accuracy of 68%, while the wet season crop was 75% accurate. The dry season rice kappa was 0.49, while the wet season was 0.60. Misclassification was commonly due to rain-fed pixels being erroneously classified as non-rice crops along with irrigated rice mistaken for rain-fed rice. When the two rice classes are combined,  the class accuracies mostly improved, as shown in Table 4. The overall accuracy improved to 76% (Dry season) and 83% (Wet season). The dry season kappa degraded slightly to 0.44, while improved for the wet season to 0.64. Rice class accuracies were high, ranging from 74% to 94%.   We compared the annual rice area at the municipal level from official statistical sources with estimates from the RS-based map ( Figure 11). The Pearson correlation coefficient was 0.75 and P < 0.001 (n = 44), indicating a strong relationship and significant correlation. RMSE was at 2,269 ha. The results show that the RS-based mapping is comparable to official government statistics at the sub-provincial scale.

Discussion
Our findings suggest that we can map irrigated and rain-fed rice with good accuracy and even better results when mapped at the rice and non-rice classes. Such maps would provide good reference data for assessing development and policy needs in the region, considering climate change and population growth issues. The design is easily transferable to similar areas with minimal reliance on ground reference data due to open data and open-source software.
One of the benefits of our approach was to use a temporal compositing method to generate seasonal statistics that served as the primary inputs. The Sentinel-2 constellation consists of two satellites with slightly different orbits; thus, the number of observations in each pixel's raw time-series will be nonuniform, and image footprints can overlap ( Figure 2). Composites such as mean, standard deviation, and quantiles can reduce these potential causes of inaccuracies. Another is that optical time series can be noisy due to atmospheric variations and observation geometries. The compositing method over a temporal window essentially eliminated noisy pixels. In some instances, statistical modeling is employed to reconstruct time-series, but in our study, where timing metrics were not vital, this is irrelevant (Hird and McDermid 2009). Based on this, mapping accuracy was aided by these statistical composites to reduce uncertainties from sensor characteristics and noise.
We used several knowledge bases to refine the supervised classification map. The presence and absence of irrigation are challenging to discern directly from imagery, particularly at large scales. Thus, ancillary information is used in addition to spectral information (Ozdogan et al. 2010). The supervised classification step yielded a rice and non-rice map. Still, with other data such as slope, elevation, and irrigation canal systems, it was possible to further separate rice into its rain-fed and irrigated ecosystems. For instance, an inspection in Google Earth of the irrigation network revealed the upper limit of elevation for irrigated rice at 50 m. Due to the need to periodically inundate paddies, the literature also suggested an upper limit of 3° slope (Brouwer et al. 1988). However, we observe a large discrepancy between official statistics and our estimates primarily due to incomplete ancillary information on irrigation networks. Only the national irrigation system (NIS) extant in the province was available, which has a service area of 27,907 ha (National Irrigation Administration 2020). However, smaller-scale systems exist as well but with numerous service areas, such as communal irrigation systems (CIS) and private irrigation systems (PIS). Their locations were unaccounted for in this study due to a lack of data. Aside from this, another reason for the lower estimates in remotely sensed data was due to the predominance of small paddy fields in the province, making it difficult for Sentinel-2 to map these accurately.
In light of these benefits, acquiring acceptable accuracy for the region is still challenging, as shown by our accuracy metrics. Overall accuracies were mainly at low to moderate levels, and even when aggregating to rice and non-rice, only a slight improvement was observed. Kappa was also modest in 0.49 (Dry season) and 0.60 (Wet season) for rain-fed/irrigated maps and not much improvement when aggregated to rice/nonrice. Despite this, our approach obtained an accuracy of up to 94% for rice and is in good agreement with official statistics at the municipal level with R 2 = 0.75 ( Figure 11). However, some of the uncertainties here lie with municipal statistics. Overall accuracy was affected by lower-class accuracies for non-rice crops. The area is intercropped with corn and various vegetables, which may confuse the classifier as they may exhibit similar temporal profiles to rice. Studies using optical imagery tend to have lower accuracies for regions like the Philippines, which are frequently cloudy and have highly fragmented rice fields, as shown by the study of Xiao et al. (2006), which yielded an R 2 of 0.6. Asilo et al. (2014) acquired an 87.2% overall accuracy (0.62 kappa) for a mainly irrigated and uniform rice area using MODIS. There were similar results of 69%-82% overall accuracy from a study on irrigated and rain-fed rice in South Asia using unsupervised ISODATA and ideal spectral signatures matching techniques with MODIS data (Gumma et al. 2016). Manjunath et al. (2015) employed the same methods and SPOT-VGT data. Still, with different knowledge bases for ecosystem types; however, there agreement at the province-level for the Philippines was low (R 2 = 0.38). Our results were able to obtain a reasonable degree of agreement at the sub-provincial (municipal) level. However, acquiring 100% accuracy is difficult due to the lack of irrigation data, cloud contamination, and small parcel sizes in Iloilo. In addition to these, interviews with the farmers revealed frequent changes to farming practices due to climate change, crop switching, or lack of capital.
We observed a higher rice area for the dry season than the wet season, contradicting official statistics. The weak performance during the dry season can be attributed mainly to confusion with non-rice crops, corroborated by the confusion matrices (Tables 3 and  Tables 4). Due to the lack of water during the dry season, the cultivation area of rain-fed rice decreases, while non-rice crops' area increases. This was not the case in the study. Since alternate crops are planted concurrently have similar durations to rice, this exacerbates confusion between the two. Additional data for more robust differentiation between rain-fed rice and non-rice crops during the dry season is needed, such as additional phenology metrics or alternative vegetation indices.
Given our results, we propose several approaches to improve accuracy. First, more variables can be used in addition to the current quantile and central tendency statistics. Other studies have used metrics such as start, end, or length of season and estimate dates of harvest and maturity (Guan et al. 2016;Jönsson and Eklundh 2004;Boschetti et al. 2009). Second, sensor fusion has proven to help combine the strengths of different satellites. This study necessitated the use of compositing to address persistent cloud contamination. Combining the rich spectral information of optical sensors with the cloud-penetrating abilities of synthetic aperture radar (SAR) platforms can potentially produce better results (Park et al. 2018;Zhang et al. 2018;Son et al. 2016).
Additionally, novel techniques using deep learning algorithms can utilize the spectral, spatial, and temporal characteristics of multi-source imagery, including optical, SAR, LiDAR, and even street view data. Several studies have shown its successful application to urban areas and crop-type delineation, potentially applicable to mapping the various rice ecosystems (Shao, Wu, and Li 2021;Shao et al. 2019;Shao, Zhang, and Wang 2017;Prins and Van Niekerk 2020). Finally, the irrigation infrastructure covered was only NIS. As geospatial inventories of these infrastructures become available, this would significantly improve estimates of irrigated and rain-fed areas. Also, climate or precipitation data can be used as an additional source of information. Due to the lack of granular data of this kind in the region, downscaled model outputs may also improve estimates.

Conclusions
Our study demonstrated the ability of a time-series technique within GEE to discriminate rice and nonrice areas in a province in the Philippines at 10 m spatial resolution. This reduces the burden of computing resources on the user, which would otherwise be challenging to scale up to more extensive regions such as in the region-wide scopes. Our method used Sentinel 2A/B, various spectral indices, and temporal composites to achieve results comparable to SARbased studies and improve upon optical-based studies. A distinctive feature of our study was the use of GCVI as the primary index for mapping rather than NDVI. The higher dynamic range of this approach enables easy detection of naturally dense vegetation such as forests and good separation of vegetation from water bodies and built-up areas. GCVI is typically used for yield or chlorophyll estimation. Still, its features allow it to be effective for mapping, mainly when temporal composites are to be used, such as annual mean or median GCVI. Another novelty in our method was to use an ensemble classification approach using the strengths of both unsupervised, supervised, and knowledge-based methods. This study could map several land cover types and rainfed and irrigated rice from this approach. We achieved a good correlation between our satellitebased estimates and official statistics at the municipal level with an R 2 = 0.75. On the provincial scale, this study yielded an 11% deviation from official reports. As this method can be observed over the same period of official harvest surveys (dry and wet seasons), this can complement those. The relative ease of access and processing of data within a cloud computing platform makes it possible for our approach to be adapted to other areas or even scaled to more extensive regions while providing time-series-based analyses.