Mapping Paddy Fields in Japan by Using a Sentinel-1 SAR Time Series Supplemented by Sentinel-2 Images on Google Earth Engine

: Paddy ﬁelds play very important environmental roles in food security, water resource management, biodiversity conservation, and climate change. Therefore, reliable broad-scale paddy ﬁeld maps are essential for understanding these issues related to rice and paddy ﬁelds. Here, we propose a novel paddy ﬁeld mapping method that uses Sentinel-1 synthetic aperture radar (SAR) time series that are robust for cloud cover, supplemented by Sentinel-2 optical images that are more reliable than SAR data for extracting irrigated paddy ﬁelds. Paddy ﬁelds were provisionally speciﬁed by using the Sentinel-1 SAR data and a conventional decision tree method. Then, an additional mask using water and vegetation indexes based on Sentinel-2 optical images was overlaid to remove non-paddy ﬁeld areas. We used the proposed method to develop a paddy ﬁeld map for Japan in 2018 with a 30 m spatial resolution. The producer’s accuracy of this map (92.4%) for non-paddy reference agricultural ﬁelds was much higher than that of a map developed by the conventional method (57.0%) using only Sentinel-1 data. Our proposed method also reproduced paddy ﬁeld areas at the prefecture scale better than existing paddy ﬁeld maps developed by a remote sensing approach.


Introduction
Rice is a staple food for billions of people especially in Asia, and paddy fields play important environmental roles by regulating water and energy budgets and supporting local biodiversity [1][2][3]. Paddy fields have also cultural and esthetic values, as several traditional ones were awarded as Globally Important Agricultural Heritage [4]. In addition, paddy fields are remarkable with respect to global warming because they are a major source of methane (CH 4 ), the second most influential greenhouse gas for global warming following CO 2 . According to the fifth assessment report of the Intergovernmental Panel on Climate Change, global emissions from paddy fields are estimated to be 33 to 40 Tg (CH 4 )/yr, or approximately 12% of total anthropogenic methane emissions [5]. In Japan, however, methane emissions from paddy fields account for more than 45% of total anthropogenic methane emissions [6]. Thus, regional methane budgets strongly depend on the paddy field distribution. High-accuracy mapping of broad-scale paddy fields is of fundamental importance for understanding these issues related to rice and paddy fields.
Satellite remote sensing is often the most appropriate approach for a comprehensive understanding of land cover. For example, the MCD12Q1 is a global land cover map based on Moderate Resolution Imaging Spectroradiometer (MODIS) data, and the Global Food Security-support Analysis Data Product 2 of 17 (GFSAD1000) is a global cropland map based on multi-satellite data (e.g., Landsat, MODIS) [7,8].
With the increasing necessity of broad-scale paddy field maps in various research fields, many efforts to produce such maps have used satellite remote sensing images, including both optical images and synthetic aperture radar (SAR) images. Before the early 2010s, most research on broad-scale paddy field mapping used remote sensing images with low to medium spatial resolutions, such as MODIS (250 to 1000 m spatial resolutions) [9][10][11][12][13][14][15][16]. Processing of multi-source remote sensing images at high spatial resolutions, such as Landsat (30 m spatial resolution), on a broad scale, is a major challenge because it requires a high computer operating capacity and large amounts of storage. Until recently few remote sensing images have been available free of charge or at a low price, broad-scale paddy field mapping has been very expensive. For this reason, various mapping techniques based on limited remote sensing resources have been developed to improve the accuracy of paddy field mapping. For example, the land surface water index (LSWI), normalized difference vegetation index (NDVI), and enhanced vegetation index (EVI), based on multispectral optical images, have been effectively used to extract irrigated paddy fields [9][10][11][12][13][14], and optical images from multiple sources have also been used to improve classification accuracy [15,16].
In recent years, many remote sensing analyses have utilized cloud-based platforms [17][18][19]. In particular, Google Earth Engine (GEE) is a popular cloud-based platform for broad-scale geospatial analyses. With GEE, high-performance computing resources for the processing of very large geospatial datasets can be accessed without the necessity to provide in-house information technology support [20]. Therefore, GEE solves many difficulties associated with broad-scale remote sensing analyses and makes it possible to use huge numbers of multisource remote sensing images with high spatial resolution. For example, Dong et al. [17], using over a thousand scenes of Landsat 8 images on GEE, have developed a paddy field map for northeastern Asia with a 30 m spatial resolution. The spatial resolution of this paddy field map is much finer than that of conventional paddy field maps developed by using low spatial resolution satellite images.
However, only a few broad-scale paddy field maps have been developed to date because some serious technical issues remain. For example, rice cultivation periods roughly overlap with rainy and cloudy seasons, and they differ among regions depending on factors such as the rice cultivar and the water management practice used. To solve these issues, several methods that use SAR images have been proposed, because SAR images can be obtained even under cloudy conditions. Among the many kinds of SAR images, Sentinel-1 images are used often because they are open source and have high spatial and temporal resolution [19,[21][22][23][24][25][26][27][28][29][30]. However, methods that use only Sentinel-1 images have low accuracy for paddy field extraction, because the backscatter of Sentinel-1 is less sensitive to vegetation and irrigated agricultural areas than multi-band optical sensors [21].
In this study, to compensate for this weakness of Sentinel-1 time series images for broad-scale paddy field mapping, we propose a novel paddy field mapping method that supplements the use of Sentinel-1 time series images with Sentinel-2 optical images to reduce misclassification of land covers. The proposed method is more reliable and robust for the extraction of irrigated areas across various land surface types. We developed two sets of paddy field maps with a 30 m spatial resolution for Japan in 2018. We developed one set by using a conventional method based on only Sentinel-1 time series images, and we developed the other set by using our proposed method based on Sentinel-1 time series images supplemented by Sentinel-2 images. Then, by comparing the results obtained by the two methods, we showed that our proposed method improved the paddy field mapping accuracy. Then, we compared the paddy field maps developed by our method with two major existing paddy field maps to demonstrate the potential of our method for accurate broad-scale paddy field mapping.
Until 2017, MAFF published typical rice cultivation dates for each prefecture, including "Transplant Start " (TS), "Transplant End" (TE), and "Harvest End" (HE) dates [33], where TS and TE are the dates on which 5% and 95%, respectively, of the total paddy field area has been transplanted, and HE is the date on which 95% of the total paddy field area has been harvested. We used the average day of the year (DOY) calculated using the TS, TE, and HE dates for the three years from 2015 to 2017 in our analysis (Figure 1).  Sentinel-1 Ground Range Detected (GRD) images were used as fundamental data for our paddy field mapping method because the C-band SAR instrument of Sentinel-1 can obtain data through cloud cover. Sentinel-1 GRD images are thus particularly appropriate for paddy field mapping in Japan, where the rainy season overlaps with the rice cultivation period. Characteristics of the satellite data used in this study are shown in Table 1. The Sentinel-1 satellites were launched by the European Space Agency (ESA) on 3 April 2014 (Sentinel-1A) and 25 April 2016 (Sentinel-1B) [34]. The SAR sensors onboard the Sentinel-1 satellites acquire C-band images at 5.4 GHz in four modes with different spatial resolutions and swath widths. We used the images acquired in Interferometric Wide Swath mode, which acquires images with a swath width of 250 km and a 10 m spatial resolution. To minimize the effects of orbit direction and look angle, for all prefectures except Okinawa prefecture, we used only images that were acquired by the Sentinel-1A satellite during descending orbits in our analysis. For Okinawa prefecture, we used images acquired by the Sentinel-1B satellite during ascending orbits. We used only VH (vertical transmit and horizontal receive) polarized images because some studies have shown that VH-polarized backscatter is more sensitive to rice growth and phenology than other polarizations [22,23]. Both Sentinel-1 satellites have a 12-day revisit cycle at the equator. Therefore, most of the study area is observed about 30 times per year. Using GEE, we accessed all Sentinel-1 GRD images (1073 images) covering the period from 1 January 2018, to 31 December 2018. Sentinel-2 Multispectral Imager (MSI) images were used in a supplemental role for the extraction of paddy fields because they are more reliable and robust for extracting irrigated areas. The Sentinel-2 satellites were launched by the ESA on 23 June 2015 (Sentinel-2A) and 7 March 2017 (Sentinel-2B) [34]. Both satellites have a 10-day revisit cycle. Therefore, most of the study area is observed more than 60 times per year. The Sentinel-2 satellites carry the Sentinel-2 MSI, which is an optical sensor with 12 bands for acquiring images with a swath width of 290 km and a spatial resolution of 10-60 m. We accessed Sentinel-2 MSI Level-1C Top of Atmosphere (TOA) reflectance images (49,694 images) acquired during the irrigated period (from 19 February 2018, to 25 September 2018) by both Sentinel-2A and Sentinel-2B satellites from GEE. Here, the irrigated period is defined as the period from TS to the date calculated by adding 30 days to TE.

Reference Agricultural Fields
Two classes of reference fields, "Paddy" and "Not paddy", were prepared for accuracy assessments; here, "Not paddy" reference fields are agricultural fields other than paddy fields.
For "Paddy" reference fields, 100 paddy fields were selected from each of 14 prefectures by using Google Maps Street View (Figure 2a), which is being increasingly used for acquiring ground truth data of various land surface types, including agricultural fields [35][36][37]. We selected the 14 prefectures from the 47 prefectures in Japan because Street View data available for the irrigation period in 2018 were limited. However, they effectively covered the whole area of Japan, allowing us to conduct data validation for the present purpose. Further validation including other prefectures remained for our forthcoming study. Irrigation management during the cultivation period of rice is distinctive from that used for other major crops. Therefore, fields that were irrigated and planted to rice could be selected by visual interpretation and designated as "Paddy" reference fields ( Figure 3). the 14 prefectures from the 47 prefectures in Japan because Street View data available for the irrigation period in 2018 were limited. However, they effectively covered the whole area of Japan, allowing us to conduct data validation for the present purpose. Further validation including other prefectures remained for our forthcoming study. Irrigation management during the cultivation period of rice is distinctive from that used for other major crops. Therefore, fields that were irrigated and planted to rice could be selected by visual interpretation and designated as "Paddy" reference fields (Figure 3).
For "Not paddy" reference fields, agricultural field polygons derived from high spatial resolution aerial photographs and provided by MAFF were used. Figure 2. (a) The 14 prefectures selected as "Paddy" reference field regions; 100 paddy fields from each of the 14 prefectures were selected as "Paddy" reference fields. (b) Six municipal districts in Hokkaido prefecture were selected as regions for "Not paddy" reference fields. All agricultural fields in these six districts were selected as "Not paddy" reference fields.

Figure 2.
(a) The 14 prefectures selected as "Paddy" reference field regions; 100 paddy fields from each of the 14 prefectures were selected as "Paddy" reference fields. (b) Six municipal districts in Hokkaido prefecture were selected as regions for "Not paddy" reference fields. All agricultural fields in these six districts were selected as "Not paddy" reference fields.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 28 Six municipal districts in Hokkaido prefecture (Figure 2b) were selected as "Not paddy" reference field regions. In most of Japan, most cultivated fields are paddy fields, but Hokkaido is atypical because broad areas are used to cultivate grains other than rice, such as wheat, corn, and soybeans. In the six selected districts, no fields were planted to rice in 2018, according to national statistics released by MAFF [32]. Therefore, all agricultural field polygons in the six districts were used as "Not paddy" reference fields.
A total of 1400 field polygons were used as "Paddy" reference fields, and 36,703 field polygons were used as "Not paddy" reference fields.

Existing Paddy Field Maps
Two different paddy field maps that had been developed by a remote sensing approach were used for accuracy assessment. The first map was produced by Takeuchi and Yasuoka [15,16] from For "Not paddy" reference fields, agricultural field polygons derived from high spatial resolution aerial photographs and provided by MAFF were used. Six municipal districts in Hokkaido prefecture ( Figure 2b) were selected as "Not paddy" reference field regions. In most of Japan, most cultivated fields are paddy fields, but Hokkaido is atypical because broad areas are used to cultivate grains other than rice, such as wheat, corn, and soybeans. In the six selected districts, no fields were planted to rice in 2018, according to national statistics released by MAFF [32]. Therefore, all agricultural field polygons in the six districts were used as "Not paddy" reference fields.
A total of 1400 field polygons were used as "Paddy" reference fields, and 36,703 field polygons were used as "Not paddy" reference fields.

Existing Paddy Field Maps
Two different paddy field maps that had been developed by a remote sensing approach were used for accuracy assessment. The first map was produced by Takeuchi and Yasuoka [15,16] from MODIS and Advanced Spacebourne Thermal Emission and Reflection Radiometer (ASTER) images obtained in the early 2000s. This map, referred to here as the "TY" map, covers Monsoon Asia (East Asia, Southeast Asia, and part of South Asia) with a 1 km spatial resolution. The second map is a land-use map produced by the Japan Aerospace Exploration Agency (JAXA) from Landsat 8 operational land imager images obtained between 2014 and 2016 [38]. This map, the "JAXA" map, covers the whole area of Japan with a 30 m spatial resolution, and one of its land-use classes is the "paddy field." Figure 4 shows a flow chart of the approach used by this study. First, two types of paddy field maps were developed, one by the "S-1" method and the other by the "S-1 & S-2" method. The S-1 method is a conventional paddy field mapping method that uses a decision tree approach and only Sentinel-1 time series data [25][26][27]. The S-1 & S-2 method is our proposed method, which also uses a decision tree approach and Sentinel-1 time series data, but it also uses an additional mask based on Sentinel-2 images. Paddy field maps for each of the 47 prefectures created by each method were merged into a single map for each method covering the whole land area of Japan. All of the processes used to develop the paddy field maps up to the final merging were performed in the GEE platform (https://earthengine.google.com/). The merging, accuracy assessment, and comparison processes were performed in ArcMap 10.5.1.   interquartile range of VH σ⁰ calculated from 2018 data in each of the 14 "Paddy" reference field regions. In all "Paddy" reference field regions, a local VH σ⁰ minimum was observed during the irrigated period. Subsequently, VH σ⁰ increased and reached a peak within a few months. Thereafter, during the latter part of the rice cultivation period, VH σ⁰ maintained a high value or decreased slightly. Previous studies have shown that VH σ⁰ decreases during the transplanting period, when the fields are flooded with irrigation water, and then increases during the ricegrowing period until around harvest time [23-28]. Therefore, rice phenology in broad areas of Japan can be observed in VH σ⁰ time series. The low VH σ⁰ values from January to March in

Preprocessing
To generate the backscatter coefficient (σ 0 ), the Sentinel-1 GRD images accessed on the GEE platform were preprocessed by applying five Sentinel-1 Toolbox corrections developed by the ESA Remote Sens. 2020, 12, 1622 7 of 17 in the following order: (1) orbit file application, (2) ground range detected border noise removal, (3) thermal noise removal, (4) radiometric calibration, and (5) terrain correction. To these preprocessed images, we applied a median filter with a 3 × 3 pixels moving window to reduce inherent speckle noise [24]. In broad-scale analyses of satellite remote sensing images, overlapping observation areas are inevitable. In the case of SAR data, overlapping areas with different incidence angles can produce noise in a time series analysis [19]. Therefore, we applied incidence angle processing to overlapping observation areas.
To the Sentinel-2 MSI Level-1C TOA reflectance images that we accessed on the GEE platform, we applied radiometric and geometric corrections, including orthorectification and spatial registration to a global reference system with sub-pixel accuracy. Optical remote sensing images are greatly influenced by cloud conditions. Therefore, pixels heavily affected by dense and cirrus clouds were removed from all Sentinel-2 MSI images by applying the QA60 Quality Assessment band. The temporal behavior of Sentinel-1 VH σ 0 over each "Paddy" reference field region was analyzed to determine the threshold values for the mask described in 2.5.3. Figure 5 shows the time series of second quartile (Q2) values and the first (Q1) to third (Q3) interquartile range of VH σ 0 calculated from 2018 data in each of the 14 "Paddy" reference field regions. In all "Paddy" reference field regions, a local VH σ 0 minimum was observed during the irrigated period. Subsequently, VH σ 0 increased and reached a peak within a few months. Thereafter, during the latter part of the rice cultivation period, VH σ 0 maintained a high value or decreased slightly. Previous studies have shown that VH σ 0 decreases during the transplanting period, when the fields are flooded with irrigation water, and then increases during the rice-growing period until around harvest time [23][24][25][26][27][28]. Therefore, rice phenology in broad areas of Japan can be observed in VH σ 0 time series. The low VH σ 0 values from January to March in Hokkaido, Aomori, Yamagata, Niigata, and Ishikawa prefectures are attributable to snow cover because these prefectures are located in northern Japan, where heavy snow falls in winter.   The first (Q1), second (Q2), and third (Q3) quartiles of VH σ 0 over "Paddy" reference fields were calculated for all pixels included within "Paddy" reference fields, with a 10 m spatial resolution.
Blue lines indicate the Q2 values, and blue shading indicates the Q1 to Q3 interval.

Masks Used for Provisional Extraction of Paddy Fields
"Forest Area (FA)" mask: A mask for forest area was made based on an updated version of the Hansen Global Forest Change map (v1.6, 2000-2018) [39]. Pixels with a forest area covering more than 30% of the pixel area were specified as forest pixels. On Hansen's map, forest area is defined by the canopy closure due to vegetation taller than 5 m. The Forest area masking was performed to minimize commission errors in paddy field mapping [17].
"Local Maximum & Minimum (LMM)" mask: The local maximum σ 0 value and the local minimum σ 0 value of the Sentinel-1 time series within 90-day moving windows were calculated for each Sentinel-1 image from the period between TS and HE. In this study, the thresholds of the local minimum and local maximum σ 0 values were determined by analyzing the temporal behavior of Sentinel-1 VH σ 0 values in the "Paddy" reference fields ( Figure 5). Figure 6 shows the local minimum Q3 VH σ 0 value during the irrigated period and the local maximum Q1 VH σ 0 value between the end of the irrigated period and HE for each "Paddy" reference field region. In all "Paddy" field reference regions, the local minimum Q3 VH σ 0 value was smaller than −20 dB. The local maximum Q1 VH σ 0 value in Hokkaido was slightly below −17 dB, whereas those in other areas were above −17 dB. Therefore, pixels with a local minimum VH σ 0 larger than −20 dB or a local maximum VH σ 0 smaller than −17 dB, were masked on each Sentinel-1 image. The VH σ 0 of paddy fields shows a very distinct minimum when the paddy fields are irrigated just after transplanting and a maximum during the heading stage of the rice plants [26].
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 28 value between the end of the irrigated period and HE for each "Paddy" reference field region. In all "Paddy" field reference regions, the local minimum Q3 VH σ⁰ value was smaller than -20 dB.
The local maximum Q1 VH σ⁰ value in Hokkaido was slightly below -17 dB, whereas those in other areas were above -17 dB. Therefore, pixels with a local minimum VH σ⁰ larger than −20 dB or a local maximum VH σ⁰ smaller than −17 dB, were masked on each Sentinel-1 image. The VH σ⁰ of paddy fields shows a very distinct minimum when the paddy fields are irrigated just after transplanting and a maximum during the heading stage of the rice plants [26].
"Local Variation (LV)" mask: Local variation values were calculated for each Sentinel-1 image by using the following equation: Local variation = Local maximum VH σ⁰ -Local minimum VH σ⁰ (1) Figure 6. The local minimum Q3 value during the irrigated period and the local maximum Q1 value between the end of the irrigated period and HE in each "Paddy" reference field region.
Local maximum VH σ⁰ and local minimum VH σ⁰ are the values calculated for the LMM mask. Pixels with a local variation of less than 5 dB were masked on each Sentinel-1 image. This threshold value, which was also determined by analyzing the temporal behavior of Sentinel-1 VH σ⁰ values in the "Paddy" reference fields, was selected to differentiate vegetation with clear seasonality from land covers with a constant low or high backscatter such as waterbodies and urban areas [26].
It should be noted that the FA mask was applied to the entire Sentinel-1 time series, whereas the LMM and LV masks were applied to each Sentinel-1 image composing the Sentinel-1 time series. The pixels of Sentinel-1 images during the irrigated period that were not masked by any of these three masks were extracted as provisional paddy field pixels. In this paper, we refer to the application of this series of three masks as the "S-1" method. We used the S-1 method to prepare a paddy field map of Japan for later comparison with a paddy field map developed by our proposed method. "Local Variation (LV)" mask: Local variation values were calculated for each Sentinel-1 image by using the following equation: local maximum VH σ 0 and local minimum VH σ 0 are the values calculated for the LMM mask. Pixels with a local variation of less than 5 dB were masked on each Sentinel-1 image. This threshold value, which was also determined by analyzing the temporal behavior of Sentinel-1 VH σ 0 values in the "Paddy" reference fields, was selected to differentiate vegetation with clear seasonality from land covers with a constant low or high backscatter such as waterbodies and urban areas [26]. It should be noted that the FA mask was applied to the entire Sentinel-1 time series, whereas the LMM and LV masks were applied to each Sentinel-1 image composing the Sentinel-1 time series. The pixels of Sentinel-1 images during the irrigated period that were not masked by any of these three masks were extracted as provisional paddy field pixels. In this paper, we refer to the application of this series of three masks as the "S-1" method. We used the S-1 method to prepare a paddy field map of Japan for later comparison with a paddy field map developed by our proposed method.

Mask Based on Sentinel-2 MSI Environmental Indexes
"Sentinel-2 indexes (S2I)" mask: The provisional paddy fields extracted by the S-1 method were reexamined by using water and vegetation indexes, LSWI, NDVI, and EVI, calculated from Sentinel-2 MSI images acquired during the irrigated period. The LSWI is sensitive to the total amount of liquid water in vegetation and its soil background [40], and the NDVI and EVI are sensitive to variations in the vegetation [41,42]. These three indexes were calculated by using Equations (2)-(4), respectively: where Blue, Red, NIR, and SWIR are the reflectance values at the top of the atmosphere of the B2 band (central wavelength: 490 nm), B4 band (665 nm), B8 band (842 nm), and B11 band (1610 nm) of the Sentinel-2 MSI. Previous studies have revealed that relationships among the LSWI, NDVI, and EVI can be used effectively to extract irrigated areas [9,17,43]. After the three S-1 masks were applied to each Sentinel-1 image, the local maximum "LSWI−NDVI" and "LSWI−EVI" values within the 10 days after the date that the Sentinel-1 image was acquired were calculated, and pixels in which the local maximum "LSWI−NDVI" and "LSWI−EVI" values were both smaller than 0 were masked.
These threshold values are based on values proposed by Dong et al. [17]. The pixels of Sentinel-1 images acquired during the irrigated period that were not masked by any of the four masks were extracted as paddy fields. We call the application of this series of four masks the "S-1 & S-2" method in this paper.

Accuracy Assessment
The accuracy of the paddy field maps was assessed by producer's accuracy with a pixel-scale resolution of 30 m. The producer's accuracy indicates the quality of the classification and is calculated by dividing the number of correctly classified reference pixels in each category by the total number of reference pixels in that category [44]. All pixels of the paddy field maps spatially included in the two types of reference field polygons were assessed for accuracy.

Accuracy Assessments in Reference Agricultural Fields
The producer's accuracy for "Paddy" reference fields was 83.6% by the S-1 method and 79.2% by the S-1 & S-2 method. For "Not paddy" reference fields, the producer's accuracy was 57.0% by the S-1 method and 92.4% by the S-1 & S-2 method. Thus, the application of the S2I mask following the three S-1 masks caused the producer's accuracy for "Paddy" reference fields to decrease by 4.4 points and that for "Not paddy" reference fields to increase by 35.4 points.
The producer's accuracies in each "Paddy" and "Not paddy" reference field region are shown in Tables 2 and 3, respectively. The producer's accuracies in "Paddy" reference field regions were 68.6-100.0% by the S-1 method and 68.1-98.3% by the S-1 & S-2 method. The producer's accuracies in "Not paddy" reference field regions were 29.6-78.4% by the S-1 method and 88.0-97.6% by the S-1 & S-2 method. Thus, the application of the S2I mask decreased the producer's accuracy in "Paddy" reference field regions by 0-16.3 points and increased the producer's accuracy in "Not paddy" reference field regions by 14.8-61.1 points.  Table 3. Producer's accuracies in each of the six "Not paddy" reference field regions in Hokkaido prefecture (30 m pixel-scale resolution). See Figure 2b for reference field region locations.

Comparison with Existing Paddy Field Maps
The paddy field maps developed by the S-1 and the S-1 & S-2 methods are compared with the JAXA and TY paddy field maps in Figures 7 and 8. The maps prepared by the S-1 and the S-1 & S-2 methods and the JAXA map have a spatial resolution of 30 m, and the TY map has a 1 km spatial resolution. The paddy field maps developed by the S-1 and the S-1 & S-2 methods have some noticeable characteristics. Compared with the JAXA and TY maps, the map produced by the S-1 method has more paddy field pixels in Hokkaido prefecture. However, the map produced by the S-1 & S-2 method has many fewer paddy field pixels in eastern Hokkaido prefecture. In general, the paddy field distribution in Japan on the maps produced by the S-1 & S-2 method is more similar to the paddy field distributions on the JAXA and TY maps. In eastern Hokkaido, rice is not grown, but non-rice grains are widely grown. Therefore, the application of the S-1 & S-2 method reduces the misclassification of non-paddy fields in that region.
We next compared the paddy field area at the prefecture scale, calculated for each paddy field map, with paddy field area data published by MAFF [32] (Figure 9). When the total paddy field area in Japan according to MAFF was set to 100%, the total paddy field area on the map produced by the S-1 and S-1 & S-2 methods was 156.7% and 112.4%, respectively, and on the JAXA and TY maps, it was 168.5% and 84.7%, respectively. The maps produced by the S-1 and the S-1 & S-2 methods greatly overestimated the paddy field area in Hokkaido prefecture (Figure 9a,b, respectively), but in other prefectures, paddy field areas did not differ significantly from the MAFF values.
To examine the overall trend of paddy field areas on each map relative to the MAFF data, we fitted straight lines to the data shown in Figure 9 by the least-squares method with the y-intercept fixed at 0. Table 4 shows the slopes of the fitted lines together with the coefficients of determination (R 2 ) and the mean square error (MSE) of the map data relative to the MAFF data. The paddy field areas calculated from the maps produced by the S-1 and S-1 & S-2 methods and the JAXA map (slopes > 1) were overestimated compared with the published MAFF areas. In contrast, the paddy field areas on the TY map (slope = 0.86) were underestimated. When Hokkaido was excluded, the areas obtained by the S-1 & S-2 method were still overestimated, but the slope was closer to 1. R 2 is a measure of how well the areas calculated from each map replicate the MAFF areas at the prefecture scale, where R 2 = 1 indicates a perfect match between them. For all of Japan, R 2 values for the S-1 and S-1 & S-2 methods were lower than those for JAXA and TY, whereas when Hokkaido was excluded, R 2 value for the S-1 & S-2 method was much higher than those for JAXA and TY. MSE is the average of the squared errors between the paddy field area in each prefecture calculated from the paddy field maps and the MAFF areas. When Hokkaido was excluded, MSE values for the S-1 & S-2 method were lower than those of JAXA and TY.   Table 4. Slopes of straight lines fitted to the relationships shown in Figure 9 with the y-intercept fixed at 0, together with the coefficients of determination (R 2 ) and the mean square errors (MSE) relative to the MAFF data.  To examine the overall trend of paddy field areas on each map relative to the MAFF data, we fitted straight lines to the data shown in Figure 9 by the least-squares method with the y-intercept fixed at 0. Table 4 shows the slopes of the fitted lines together with the coefficients of determination (R²) and the mean square error (MSE) of the map data relative to the MAFF data. The paddy field areas calculated from the maps produced by the S-1 and S-1 & S-2 methods and the JAXA map (slopes > 1) were overestimated compared with the published MAFF areas. In contrast, the paddy field areas on the TY map (slope = 0.86) were underestimated. When Hokkaido was excluded, the areas obtained by the S-1 & S-2 method were still overestimated, but the slope was closer to 1. R² is a measure of how well the areas calculated from each map replicate the MAFF areas at the prefecture scale, where R² = 1 indicates a perfect match between them. For all of Japan, R² values for the S-1 and S-1 & S-2 methods were lower than those for JAXA and TY, whereas when Hokkaido was excluded, R² value for the S-1 & S-2 method was   Table 4. Slopes of straight lines fitted to the relationships shown in Figure 9 with the y-intercept

Discussion
We proposed a new high-resolution paddy field mapping method that uses Sentinel-1 SAR time series images supplemented by Sentinel-2 optical images. Although some paddy field mapping methods combining SAR and optical images have been proposed previously [19,21,22], those methods used optical images as fundamental data for paddy field extraction; therefore, they could be seriously affected by gaps and biases in the optical images due to cloud cover. In contrast, in our proposed method, the Sentinel-2 images are used only as ancillary data to produce masks based on water and vegetation indexes. As a result, gaps in the Sentinel-2 images did not adversely affect the paddy field mapping accuracy. In this regard, our method is expected to have an advantage for detecting paddy fields at a broad scale.
The producer's accuracy for "Paddy" reference fields identified by the S-1 & S-2 method was lower than that of fields identified by the S-1 method; this result is reasonable because the additional mask applied to extract paddy fields in the S-1 & S-2 method makes the identification criteria more strict. In contrast, the producer's accuracy for "Not paddy" reference fields identified by the S-1 & S-2 method was higher compared with the S-1 method. Therefore, to grasp the potential of the S-1 & S-2 method for paddy field mapping, a comprehensive evaluation is needed to account for the decrease of the producer's accuracy in "Paddy" reference fields and the increase in "Not paddy" reference fields. The map developed by the S-1 & S-2 method has a spatial resolution of 30 m, which is much finer than previous maps based on satellite images with low to medium spatial resolutions. However, this resolution might be still insufficient for validation of the accuracy at the pixel scale for individual agricultural fields in Japan. In 2009, 38% of paddy fields in Japan were smaller than 0.003 km 2 , and 92% of paddy fields were smaller than 0.01 km 2 , according to MAFF [45]. Therefore, we believe that the producer's accuracy for "Paddy" reference fields (79.2%) identified by the S-1 & S-2 method is still acceptable, and the decrease of 4.4 points in the accuracy for "Paddy" reference fields compared to the S-1 method does not indicate a serious reduction in paddy field extraction accuracy at a broad scale. In contrast, the accuracy for "Not paddy" reference fields of the S-1 & S-2 method was greatly improved, by 35.4 points, compared with that of the S-1 method. This large increase in accuracy justifies the use of the S-1 & S-2 method despite the decrease in accuracy for "Paddy" reference fields. Therefore, we consider the S-1 & S-2 method to be an effective paddy field mapping approach because it reduces the misclassification of non-paddy agricultural fields as paddy fields.
Remarkably, the effectiveness with which the S-1 & S-2 method reduced misclassification differed among regions (Table 4). In particular, in the Koshimizu district, Hokkaido prefecture, 12% of "Not paddy" reference fields were misclassified. The inter-regional variation in accuracy is primarily attributable to gaps in Sentinel-2 images caused by cloud cover. Although unlike previous paddy field mapping methods that use optical images as fundamental data, the S-1 & S-2 method does not generally misclassify or produce no-data pixels because of gaps in the Sentinel-2 images, if there are no cloud-free Sentinel-2 images during the irrigated period, the S2I mask is ineffective. In this situation, therefore, gaps in the Sentinel-2 images can reduce the degree to which the classification accuracy is improved by the S-1 & S-2 method. It is difficult to completely resolve this issue because the effects of cloud cover cannot be controlled, but it should be possible to reduce the uncertainty of the classification by using a larger number of cloud-free optical image scenes. Fortunately, the ESA has planned the launches of Sentinel-2C and Sentinel-2D, and these additional satellites will make it possible to use more Sentinel-2 images in the future. This future increase in the number of Sentinel-2 satellites will reduce the uncertainty of the S-1 & S-2 method by improving its temporal resolution.
This study developed a paddy field map covering a wide area of Japan, except for Hokkaido prefecture, that reproduced well data from MAFF on paddy field areas. The maps produced by the S-1 and S-1 & S-2 methods, however, seriously overestimated the paddy field area in Hokkaido prefecture ( Figure 9). This prefecture is in a cooler climate and has more types of agricultural fields covering a larger area (11,450 km 2 ) than the other prefectures; in 2018, the cultivated paddy field area accounted for only 9.3% of the total agricultural field area in Hokkaido prefecture [32]. In contrast, in the other Japanese prefectures, the cultivated paddy field area accounted for 45.4% of the total agricultural field area in 2018 [32]. The maps by JAXA and TY used multi-band optical sensors that are more robust for extracting vegetation type and irrigated agricultural areas in the Hokkaido prefecture than the SAR sensor. One approach to solving the overestimation in Hokkaido prefecture is, as described above, to improve the accuracy of the S2I mask for extracting irrigated areas from satellite images by using a larger number of optical image scenes. Another approach is to refine the threshold values of the LMM and LV masks to improve the accuracy with which paddy fields are extracted from the Sentinel-1 time series data. In this study, we used the same threshold values for these masks to extract paddy fields everywhere in Japan, but the temporal behavior of Sentinel-1 VH σ 0 differed greatly among the regions (Figures 5 and 6). Therefore, to further improve the accuracy of paddy field mapping in different regions, appropriate threshold values should be selected for each region instead of using uniform country-scale values. For Japan excluding Hokkaido prefecture, however, the S-1 & S-2 method satisfactorily reproduced MAFF paddy field areas at both country scale and prefecture scale. Therefore, we are convinced that our novel S-1 & S-2 method has a high potential for paddy field mapping in most of Japan, except for the Hokkaido prefecture.
Our efforts in this study are a first step toward providing more accurate paddy field maps wherever paddy rice is grown. Because rice cultivation method and cultivation environments vary greatly from region to region, it is important to confirm whether our method is applicable to mapping paddy fields in other parts of Asia and the world. We examined the new mapping method only in Japan and found that the new method worked well in most of the study area but overestimated in Hokkaido prefecture. This implies that the method can be improved by applying additional masks and thresholds such as a more detailed classification of cultivated crops and land cover types, and by conducting validations over broader and more diverse areas. Here, by using Google Maps Street View, it was relatively easy to select "Paddy" reference fields without on-site visits, but Street View recordings are not always available in every targeted region and time. Furthermore, many paddy fields are likely distributed in areas where few Street View recordings have been made. Therefore, in such regions, it will be necessary to use different approaches for identifying reference fields. Nevertheless, in the future, it would be desirable to expand the target region to regions outside of Japan and to verify the effectiveness of our method after overcoming the issues described above.

Conclusions
Through developing paddy field maps for Japan, this study demonstrated the potential of our novel paddy field mapping method using the Sentinel-1 time series supplemented by Sentinel-2 optical images in comparison with the conventional paddy field method using only Sentinel-1 time series data. Furthermore, by comparing our maps with existing paddy field maps and MAFF data on paddy field areas, we evaluated the reproducibility of our paddy field map. Our method, which uses an additional mask based on Sentinel-2 indexes obtained during the irrigated period, shows great potential for reducing the misclassification of non-paddy agricultural fields as paddy fields. The producer's accuracy for non-paddy reference agricultural fields increased from 57.0% with the conventional method to 92.4% with our proposed method. The paddy field map developed by our proposed method also reproduced well paddy field area data at the prefecture scale more reliably than existing paddy field maps of Japan, except in Hokkaido prefecture where our method systematically overestimated paddy field area. Our advanced method would be extended by covering other land cover types such as cropland, forest, and urban areas, and therefore methodological accuracy can be improved by conducting validations over a wide variety of areas not only in Japan but also in other countries.
Author Contributions: Conceptualization and data collection, S.I. and A.I.; methodology, analysis, validation, and paper writing, S.I.; paper revision, A.I. and C.Y. All authors have read and agreed to the published version of the manuscript.