*2.1. Sites along with the Collection of Field Data*

We prepared ground observation datasets obtained at six sites (A–E) located in six different districts: Site A, in Thot Not, Can Tho (10◦10 N, 105◦33 E); Site B, in Chau Thanh (10◦16 N, 105◦08 E); Site C, in Cho Moi (10◦25 N, 105◦27 E); Site D, in Thoai Son (10◦16 N, 105◦08 E); and Site E, in Tri Ton, An Giang (10◦23 N, 105◦06 E) [7,8,10,21–27] Figure S1. The soils at sites A–C are classified as silty clay fluvisol (a type of alluvial soil; [17]), while the soils at sites D and E are classified as sulfuric humaquepts (a type of alluvial soil [17]).

In Can Tho and An Giang, 50 farmers' rice paddies (30 in site A, five each in sites B–E) were chosen as regions of interest (ROIs). At the center of each ROI, field water level data were collected for the supervised classification of the satellite remote sensing data with a water level gauge (daily, 10:00 AM–12:00 PM at site A) or with a HOBO CO-U20L-04 water level data logger (Onset Computer Corporation, United States; collected every 4 h at sites B–E). At the same time, we collected information about the history of field operations (e.g., fertilization and land preparation/sowing/harvesting dates) at each ROI throughout the

observation period. The numbers of ROIs in the inundated/non-inundated rice paddies are described in the cited literature [10].

**Figure 1.** Flowchart of the methodology of this study outlining the data, processing and analysis steps.

*2.2. CyGNSS GNSS-R Datasets and Their Preprocessing Methods*

All CyGNSS [10] Lv.1 Version 3.0 data observed from the first observation (August 2018) until December 2021 were downloaded from https://podaac.jpl.nasa.gov/ dataset/CYGNSS\_L1\_V3.0 (accessed on 15 January 2022). The reflectivity (Γ) data were calibrated [28] using Equation (1):

$$\Gamma(\theta) = \frac{(4\pi)^2 (P\_{DDM} - N) \left(R\_r + R\_t\right)^2}{\lambda^2 G\_r G\_t P\_t} \tag{1}$$

where *PDDM* is the maximum value of the analog power in the delay/Doppler maps (DDM), *N* is the noise floor related to the *DDM*, *Rr* is the receiver–specular point (SP) distance, *Rt* is the transmitter–SP distance, *λ* is the wavelength, *θ* is the incidence angle, *Gr* is the receiver antenna gain in the direction of the SP and *GtPt* is the transmitter equivalent isotropically radiated power (EIRP). The noise floor is computed as the mean value of the DDM subset, where the signal is absent (located above the characteristic horseshoe shape of the DDMs). The effect of the scattering area with the highest analog power in the DDM maps was used as the size of the specular point. Since the CyGNSS GPS signal integration time is fixed at 1 s, the footprint shape was inversely computed using the integration time, the velocity of the SP and the effective scattering area.

Our precision index model's design was inspired by the spatial localization technique of common data assimilation methods such as the local ensemble transform Kalman filter or local particle filters [10]. The precision index (PI) was calibrated using the following equation in the grid covered by the SP footprint, as shown in Equation (2):

$$PI = \frac{\cos(\theta) \times GS}{sqrt(ESA) \times \exp(3.0 \times \left(DistSP/SemiDSP\right)^2)} \approx \text{ObsEr} \tag{2}$$

where *θ* is the incidence angle, *GS* is the grid spacing, *ESA* is the effective scattering area, *DistSP* is the distance from the center of the specular point, *SemiDSP* is the semidiameter of the ellipsoidal-shaped specular point, and *ObsEr* is the observation error (Figure 2).

**Figure 2.** Illustration of the precision index. The light blue tiles are rasterization grid cells. The yellow circular area is the effective scattering area. The red tile in the yellow/green circle effective scattering area is the corresponding grid. The green, blue and red arrows are equivalent to the *GS*, *DistSP* and *SemiDSP* terms in Equation (2).

Each specular point in the CyGNSS data format contains analog power in a 17 × 11 array of DDM bins [17 rows for Delay with a 0.25-chip resolution, 11 columns for Doppler with a 500-Hz resolution]. We also analyzed the analog power in the DDM by regarding the power as a probability density of a 3-dimensional histogram (representing skewness and kurtosis) as described in Equation (3) and Figure 3 by targeting 5 × 5 arrays surrounding the element containing the maximum analog power over the DDM at each specular point. If the kurtosis value was greater than 0.01, the precision index (*P*) zeroed out before its use to omit the noise derived from the specular effects over highly rough land surfaces.

$$\begin{aligned} SKW\\_a &= \sum\_{i=-2}^{2} \sum\_{i=-2}^{2} \left( \frac{(p\_{ij} - p\_{00}) / (TP \times 25)}{1/3 + \sqrt{2} [(DPL\\_P\_{ij} - DPL\\_P\_{00})]} \right)^3 \\ SKW\\_b &= \sum\_{i=-2}^{2} \sum\_{i=-2}^{2} \left( \frac{(p\_{ij} - p\_{00}) / (TP \times 25)}{1/3 + \sqrt{2} [- (DPL\\_P\_{ij} - DPL\\_P\_{00}) + (DLY\\_P\_{ij} - DLY\\_P\_{00})]} \right)^3 \\ SKW\\_c &= \sum\_{i=-2}^{2} \sum\_{i=-2}^{2} \left( \frac{(p\_{ij} - p\_{00}) / (TP \times 25)}{1/3 + \sqrt{2} [(DPL\\_P\_{ij} - DPL\\_P\_{00}) + (DLY\\_P\_{ij} - DLY\\_P\_{00})]} \right)^3 \\ DDM\\_SkSW &= MAX \{ abs(SKW\\_a), abs(SKW\\_b), abs(SKW\\_c) \} \\ DDM\\_KTS &= \sum\_{i=-2}^{2} \sum\_{i=-2}^{2} \left( \frac{(p\_{ij} - p\_{00}) / (TP \times 25)}{1/3 + \sqrt{2} \left[ (DPL\\_P\_{ij} - DPL\\_P\_{00})^2 + \left( Dly\\_P - DPL\\_P\\_P \right)^2 \right]} \right)^4 \end{aligned} \tag{3}$$

where *SKW* is skewness, KTS is kurtosis, *P*<sup>00</sup> is the maximum analog power value (W) on the DDM arrays of each specular point ["00" indicates the index of the element containing the maximum analog power among all arrays in the *DDM*; i.e., *P*<sup>00</sup> in Equation (3) is equivalent to *PDDM* in Equation (1)], *i* and *j* are array indexes over the *DDM* surrounding the maximum analog power element (*i* is the Delay row index and *j* is the Doppler column index), *TP* indicates the sum of the power analog values of all arrays in the *DDM*, and *DPL* and *DLY* indicate the doppler–delay index of the target element (i.e., *Pij* or *P*00) over the *DDM* array.

**Figure 3.** Illustration of DDM 3D statistics (skewness and kurtosis). The calibration is conducted by assuming that the analog power is equivalent to the probability density of the DDM 3D histogram (Delay, Doppler and Analogue power).

After calibrating the reflectivity (Γ) and PI on a latitude/longitude map (with a 500 m resolution and a daily temporal resolution) (Lv. 2, Figure 4), the data were applied to a temporal Kalman smoother on each 15-day cycle (temporal localization scale: 14 days. 1σ = 5 days) to obtain the Lv. 3 product (Figure 5) for the subsequent spatiotemporal analysis. The Lv. 2 data were also applied for the temporal analysis (with a slight modification to the change detection algorithm described in [29]) just after being applied in the Γ(*θ*)-normalization task with Equation (4) and in a 30-day moving average filter; then, the results were compared with the ALOS-2/PALSAR-2 products reported in [10]. To generate the Lv. 3 products, this study simply used a linear Kalman filter (i.e., the time-evolution of the model was assumed to be negligible). Γ-reflectivity was treated as both the states and measurements.

<sup>Γ</sup>*normalized* = <sup>Γ</sup>−Γmin <sup>Γ</sup>max−Γmin referring to a paper [29]

$$\left(\Gamma(\theta)\_{normalized} = \frac{\Gamma(\theta) - \Gamma(\theta)\_{\min}}{\Gamma(\theta)\_{\max} - \Gamma(\theta)\_{\min}}\right) \tag{4}$$

where Γ(*θ*)max/min is the temporal maximum/minimum value of the corresponding incidence angle bin. Due to the data quantity limitations of specular points obtained during the 2018–2022 period, we calibrated the incidence angle bins prepared for every 5◦ interval (i.e., 0–5◦, 5–10◦, 10–15◦, 15–20◦, 20–25◦, 25–30◦, 30–35◦, 35–40◦, 40–45◦, 45–50◦, 50–55◦, 55–60◦, 60–65◦ and 65–70◦ bins) in each grid.

**Figure 4.** A sample of the Lv. 2 daily product (**a**) Γ(*θ*)normalized; (**b**) incidence angle; (**c**) precision index clipped at the Mekong Delta ((**d**) Optical image) on 4 January 2021.

**Figure 5.** A sample of the Lv. 3 15−day−cycle Kalman smoother product based on the precision index [Γ(dB) without or with applying the precision index (**a**,**b**), zeroed out based on the kurtosis threshold and DDM 3D statistics such as skewness (**c**) and kurtosis (**d**)] clipped at Mekong Delta on 1 August 2018.

After calibrating the reflectivity (Γ) and PI on a latitude/longitude map (with a 500 m resolution and a daily temporal resolution) (Lv. 2, Figure 4), the data were applied to a temporal Kalman smoother on each 15-day cycle (temporal localization scale: 14 days. 1σ = 5 days) to obtain the Lv. 3 product (Figure 5) for the subsequent spatiotemporal analysis. The Lv. 2 data were also applied for the temporal analysis (with a slight modification to the change detection algorithm described in [29]) just after being applied in the Γ(*θ*)-normalization task with Equation (4) and in a 30-day moving average filter; then, the results were compared with the ALOS-2/PALSAR-2 products reported in [10].
