*2.2. Satellite Imagery Datasets*

#### 2.2.1. PlanetScope CubeSat Surface-Reflectance Ortho Tiles

We used the OFRs' shapefile to search for and clip Level 3A surface-reflectance imagery available through Planet Orders API. The PlanetScope surface-reflectance ortho tiles use a fixed UTM grid system in 25 km by 25 km tiles with 1 km overlap [1]. We filtered out all images with more than 10% cloud using an image-based cloud-cover filter—this cloud-cover filter threshold allowed us to download mostly cloud-free images; however, because it is an image-based filter rather than an OFR or area-of-interest-based cloud filter, some useful observations (i.e., when the OFR is not covered with clouds but the image is filtered out) were not downloaded, decreasing the total number of observations per OFR. In addition, to deal with potential cloud-obscuration outliers, we used the PlanetScope

UDM2 to filter out all image clips that contained more than 5% unusable pixels (i.e., pixels covered by clouds, cloud shadow, with light and heavy haze).

The PlanetScope ortho tiles were resampled to 3 m and projected using the WGS84 datum. The ortho tiles were radiometrically, sensor, and geometrically corrected and aligned to a cartographic map projection. These images were atmospherically corrected using the 6S radiative transfer model with ancillary data from MODIS [1,38,39], and the positional accuracy has been reported to be smaller than 10 m [1].

#### 2.2.2. Normalized Surface-Reflectance Basemap

We processed daily Basemap images corrected to surface reflectance using PlanetScope scenes, and a "best scene on top" algorithm [16,18] that selects the highest quality imagery from the PlanetScope catalog. This algorithm ranks the PlanetScope scenes based on their quality by assessing the scenes' acutance (i.e., sharpness), the fraction of cloud cover, cloud shadow, haze, and presence of unusable pixels (e.g., no data). Briefly, this algorithm is based on a linear regression model approach that uses the clear pixels from the best-ranked scenes; we selected the best scenes first, then progressed successively until the images were filled or no scenes remained [18]. To obtain Basemap at a daily cadence, we employed a 30-day rolling window that may use PlanetScope scenes collected up to 15 days before or after the target date; however, if no usable pixels (i.e., cloud-free) are available in this time range, the image will contain no data. We did not observe any Basemap image with no-data in our study period. The rolling window approach weights on the image recency, for instance, a slightly hazy scene (e.g., ~<10% hazy pixels) on the day of the Basemap image, will score higher than a very clear scene (i.e., no haze) from a few days before/after. In addition, due to the daily cadence, there may be Basemap images with the same PlanetScope scene composition, which leads to repetitive information when using the Basemap images to monitor OFRs.

Basemap images were generated employing a two-step process: normalization and seamline removal. Normalization aims to radiometrically calibrate the Basemap images and to minimize the scene-to-scene variability when mosaicking PlanetScope scenes. For this step, the Framework for Operational Radiometric Correction for Environmental Monitoring (FORCE) [40] was used to generate a combined Landsat 8 and Sentinel-2 surface-reflectance product to be used as the "gold" radiometric reference during normalization. FORCE infers surface reflectance from Landsat 8 and Sentinel-2 by employing the 5S (simulation of the satellite signal in the solar spectrum) approach [41]. The aerosol optical depth is estimated using a dark-object-based approach where in water vapor content is derived from Landsat 8 (obtained from MODIS database) and Sentinel-2 (estimated on a pixelspecific basis) imagery. In addition, clouds and shadows are detected using a modified version of Fmask [42] for Sentinel-2 images [43] (see [16,17,44] for further details). An assessment of the FORCE atmospheric correction was performed as part of the atmospheric correction inter-comparison exercise [45], and the FORCE implementation uses the Landsat 8 and Sentinel-2 imagery mapped onto a common UTM grid to produce 30 m spatialresolution imagery. Seamline removal enhances the visual appearance of the Basemap image edges. In this step, each PlanetScope scene used in the Basemap mosaic is set to match its neighbor—pixel values near a scene boundary change more than values away from the boundary; however, the pixel values are not modified. Specifically, we first calculated the Basemap mosaic pixel values gradient, then set the gradient values between 1 and 0 (scene boundary) and fixed the original pixel values along the Basemap mosaic edge. This process was applied independently for each band; therefore, it may alter band ratios near scene edges—this is most apparent when scenes do not match locally, for instance, for unmasked clouds. Lastly, the seamline removal may introduce artifacts (e.g., straight lines, distortions) at the Basemap mosaic boundary, which is most frequent over water when normalization cannot fully correct for differences between scenes due to waves and sun glint.

#### 2.2.3. Planet Fusion Surface Reflectance

We processed Planet Fusion images using an algorithm based on the CubeSat-enabled spatiotemporal enhancement method [8], which enhances, inter-calibrates, and fuses satellite imagery from multiple sensors. Planet Fusion has unique features, including (1) precise co-registration and sub-pixel fine alignment for different image sources, (2) PlanetScope scenes with near-nadir field of view, resulting in minimal bidirectional reflectance distribution function (BRDF) variation effects, and (3) pixel traceability to identify imagery sources and to assess the confidence of daily gap-filled images.

To generate Planet Fusion surface-reflectance images, we used the same approach described for Basemap (i.e., FORCE [40]), with top-of-atmosphere (TOA) PlanetScope scenes (3 m), Sentinel-2 TOA reflectance (10–20 m), Landsat 8 TOA reflectance (30 m), and daily tile-based MODIS or VIIRS normalized to a nadir-view direction and local-solar-noon surface reflectance. The Planet Fusion algorithm uses MODIS MCD43A4 surface-reflectance product in seven spectral bands that are corrected for reflectance anisotropy using a semiempirical BRDF [46], which utilizes the best observations collected over a 16-day period centered on the day of interest. In addition, VIIRS products (VNP43IA4 and VNP43MA4) are used as a backup to ensure continuity if MODIS data is not available.

The Planet Fusion algorithm guarantees spatially complete and temporally continuous images by gap-filling radiometric data (i.e., synthetic pixel values). The gap-filling process uses both spatial (i.e., neighboring and class-specific pixel information) and temporal interpolation techniques to estimate the pixel values. In general, uncertainty will vary based on Earth's surface characteristics (e.g., vegetation dynamic changes), and it will be higher for longer daily interval gaps. Planet Fusion images are accompanied by a quality-assurance product, which is a thematic raster layer using the same spatial grid (i.e., UTM grid system in 24 km by 24 km tiles) as the corresponding Planet Fusion spectral data [17]. We used the quality-assurance product to assess the percentage of synthetic (i.e., gap-filled) versus observation data (PlanetScope and Sentinel-2) used to generate the pixel value. The observation data can be a combination of PlanetScope and Sentinel-2. A value of 1 indicates no gap-filling, whereas a value of 100 indicates an entirely gap-filled pixel value. Specifically for our study case, when clipping Planet Fusion images using the OFR boundaries, the clips can have real pixels, synthetic pixels, or a combination of both. Additionally, there are known issues associated with the gap-filling process used by the Planet Fusion algorithm, including false cloud or cloud-shadow detection and image artifacts (e.g., strips, distortions). These issues are most common during prolonged cloudiness and in study regions with significant terrain shadowing.

#### *2.3. Data Analysis*

To classify the OFR surface area from PlanetScope, Basemap, and Planet Fusion, we clipped all available images using each OFR shapefile buffered to 100 m. Then, we calculated the normalized difference water index (NDWI) using the green and NIR bands [47], and we applied an adaptive Otsu threshold [48] for each image in the time series to separate water from non-water pixels. The Otsu threshold is a well-known algorithm used to classify surface water of inland water bodies [2,4,49–53]. In addition, the Otsu threshold optimizes the separability of pixel values is contingent on the bimodal distribution of the pixel values (i.e., water and non-water pixels), which was ensured by clipping the satellite imagery using each OFR shapefile. After calculating the Otsu threshold and separating water from non-water pixels, we clipped the images one more time using the OFRs' shapefiles buffered at 20 m. This last step was done to minimize the impact (i.e., inflating surface area) of adjacent water bodies when estimating OFR surface area. All surface area image classification was done in Google Earth Engine [54].

PlanetScope has a near-daily revisit time; however, the number of usable satellite images varies throughout the year due to the presence of clouds and sensor-related issues. To assess the number of different PlanetScope observations, we first counted the total number of observations for each OFR and for each month; then, we plotted the monthly

distribution of this number, including all OFRs (i.e., one boxplot for each month of the year that represents the variability in the number of monthly observations according to different OFRs). In addition, we evaluated the total number of different observations for each OFR along the year (i.e., a histogram that represents the distribution of the total number of observations for each OFR). A similar approach was used to assess the number of different Basemap observations and to count the number of Planet Fusion observations that were real, mixed (i.e., including real and synthetic pixels), and synthetic. Basemap images with the same PlanetScope scene composition were counted only once.

For the different OFR surface area size classes (Figure 1), we assessed the uncertainties in the Basemap and Planet Fusion images by pairwise comparing them with PlanetScope and calculating the percent error (Equation (1)) monthly distribution and the monthly mean absolute percent error (MAPE; Equation (2)). In addition, for Planet Fusion, we divided the pairwise comparisons between real, mixed and synthetic surface area observations. We illustrated the surface area time series derived from PlanetScope, Basemap, and Planet Fusion for six OFRs of different sizes (Table 2). These OFRs were chosen to demonstrate the surface area time series variability from the different images and for OFRs located under different environmental conditions (e.g., presence of vegetation inside the OFR, close to adjacent water bodies, a multi-part OFR). In addition, we overlaid the OFRs' shapefile on high-resolution Google Maps satellite imagery to show the environmental conditions where the OFRs are located.

$$\text{Percentage}(\%) = ((y\_i - x\_i) / \chi\_i) \,\*\, 100\tag{1}$$

$$\text{Mean absolute percent error} (\%) = \frac{1}{n} \Sigma \left| \frac{y\_i - x\_i}{x\_i} \right| \ast 100 \tag{2}$$

where *xi* is the SkySat or PlanetScope surface area and *yi* is the Basemap or Planet Fusion surface area.

**OFR id OFR Size (ha)** A 13.62 B 22.60 C 9.83 D 12.73 E 13.82 F 29.72

**Table 2.** Selected OFRs to illustrate PlanetScope, Basemap, and Planet Fusion surface area time series and their size according to the OFR dataset.
