Monitoring Small Water Bodies Using High Spatial and Temporal Resolution Analysis Ready Datasets

: Basemap and Planet Fusion—derived from PlanetScope imagery—represent the next generation of analysis ready datasets that minimize the effects of the presence of clouds. These datasets have high spatial (3 m) and temporal (daily) resolution, which provides an unprecedented opportunity to improve the monitoring of on-farm reservoirs (OFRs)—small water bodies that store freshwater and play important role in surface hydrology and global irrigation activities. In this study, we assessed the usefulness of both datasets to monitor sub-weekly surface area changes of 340 OFRs in eastern Arkansas, USA, and we evaluated the datasets main differences when used to monitor OFRs. When comparing the OFRs surface area derived from Basemap and Planet Fusion to an independent validation dataset, both datasets had high agreement (r 2 ≥ 0.87), and small uncertainties, with a mean absolute percent error (MAPE) between 7.05% and 10.08%. Pairwise surface area comparisons between the two datasets and the PlanetScope imagery showed that 61% of the OFRs had r 2 ≥ 0.55, and 70% of the OFRs had MAPE <5%. In general, both datasets can be employed to monitor OFRs sub-weekly surface area changes, and Basemap had higher surface area variability and was more susceptible to the presence of cloud shadows and haze when compared to Planet Fusion, which had a smoother time series with less variability and fewer abrupt changes throughout the year. The uncertainties in surface area classiﬁcation decreased as the OFRs increased in size. In addition, the surface area time series can have high variability, depending on the OFR environmental conditions (e.g., presence of vegetation inside the OFR). Our ﬁndings suggest that both datasets can be used to monitor OFRs sub-weekly, seasonal, and inter-annual surface area changes; therefore, these datasets can help improve freshwater management by allowing better assessment and management of the OFRs.


Introduction
Planet Labs currently operates more than 200 PlanetScope satellites in sun-synchronous orbits and frequently launches new satellites that are designed to have a short operational lifetime (<4 years). The PlanetScope satellite constellation enables near-daily monitoring with multi-spectral imagery at high spatial resolution (3 m) [1]. PlanetScope imagery has been applied to a variety of studies to monitor phenomena that require both high spatial and temporal resolution, for instance, to monitor small water bodies [2][3][4], estimate methane emissions from forested wetlands [5], assess river-ice and water velocity [6], improve temporal resolution (e.g., PlanetScope [3 m] and Sentinel-2 [10 m]) when compared to Landsat. However, a multi-sensor approach requires processing of satellite imagery of different spatial resolution from multiple platforms, which can be time-consuming and a limiting factor if it is necessary to process, download, and move the satellite imagery across multiple platforms [33]. In this study, we propose a novel use of the analysis ready datasets Basemap and Planet Fusion, and we aim (1) to assess the usefulness of both datasets to monitor OFRs sub-weekly surface area changes and (2) to compare the two datasets and describe their differences when used to monitor OFRs.

Study Region
Eastern Arkansas is one of the largest irrigated regions in the USA that has seen a rapid increase in the number of OFRs during the last 40 years [35][36][37]. The region has a humid subtropical climate with an average annual precipitation of 1300 mm, mostly distributed between March and May and November and January [23]. Recent studies [35,36] mapped the spatial distribution of 340 OFRs with surface area <30 ha and distributed across three sub-watersheds in the study region ( Figure 1). The OFR dataset was manually mapped using the high-resolution (1 m) National Agriculture Imagery Program archive in combination with 2015 Google Earth satellite imagery. The authors of the OFR dataset used Google Earth Explorer to sharpen the image details when zooming in and to provide a validation for features appearing indistinct or pixelated in the 1-m mosaic imagery [35]. We assigned the OFRs to three size classes (0.1-5 ha, 5-10 ha, and 10-30 ha) based on the surface area mapped in the OFR dataset. These classes were used to support the surface area monitoring analyses when accounting for different OFR sizes ( Figure 1). Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 21 uncertainties ~ 20%). Aiming to overcome these limitations, other studies [4,33,34] have applied a multi-sensor satellite imagery approach, including sensors of higher spatial and temporal resolution (e.g., PlanetScope [3 m] and Sentinel-2 [10 m]) when compared to Landsat. However, a multi-sensor approach requires processing of satellite imagery of different spatial resolution from multiple platforms, which can be time-consuming and a limiting factor if it is necessary to process, download, and move the satellite imagery across multiple platforms [33]. In this study, we propose a novel use of the analysis ready datasets Basemap and Planet Fusion, and we aim (1) to assess the usefulness of both datasets to monitor OFRs sub-weekly surface area changes and (2) to compare the two datasets and describe their differences when used to monitor OFRs.

Study Region
Eastern Arkansas is one of the largest irrigated regions in the USA that has seen a rapid increase in the number of OFRs during the last 40 years [35][36][37]. The region has a humid subtropical climate with an average annual precipitation of 1300 mm, mostly distributed between March and May and November and January [23]. Recent studies [35,36] mapped the spatial distribution of 340 OFRs with surface area <30 ha and distributed across three sub-watersheds in the study region ( Figure 1). The OFR dataset was manually mapped using the high-resolution (1 m) National Agriculture Imagery Program archive in combination with 2015 Google Earth satellite imagery. The authors of the OFR dataset used Google Earth Explorer to sharpen the image details when zooming in and to provide a validation for features appearing indistinct or pixelated in the 1-m mosaic imagery [35]. We assigned the OFRs to three size classes (0.1-5 ha, 5-10 ha, and 10-30 ha) based on the surface area mapped in the OFR dataset. These classes were used to support the surface area monitoring analyses when accounting for different OFR sizes ( Figure 1).  We downloaded PlanetScope images and processed daily Basemap and Planet Fusion images between July 2020 and July 2021. This time frame was chosen based on the imagery availability to generate both analysis ready datasets. The images spatial resolution and band-wavelength ranges are presented in Table 1. In addition, the general workflow used to assess the OFRs' surface area time series is provided in Figure 2. We downloaded PlanetScope images and processed daily Basemap and Planet Fusion images between July 2020 and July 2021. This time frame was chosen based on the imagery availability to generate both analysis ready datasets. The images spatial resolution and band-wavelength ranges are presented in Table 1. In addition, the general workflow used to assess the OFRs' surface area time series is provided in Figure 2. 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).   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 Remote Sens. 2021, 13, 5176 5 of 20 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].

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.

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.

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][50][51][52][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.
where x i is the SkySat or PlanetScope surface area and y i is the Basemap or Planet Fusion surface area.

Validation Scheme
To validate the surface area classification using the Otsu thresholding approach, we downloaded five orthorectified and multispectral SkySat images [17] (Blue: 0.450-0.515 µm, Green: 0.515-0.595 µm, Red: 0.605-0.695 µm, and NIR: 0.740-0.900 µm) at sub-meter (0.66-0.73 m) spatial resolution (Table 3). For each image, we overlaid the OFR geometry and manually delineated the OFR surface area, which resulted in 144 validation surface areas from 71 different OFRs for multiple observations in time. Then, we conducted a pairwise comparison of the validation surface area with the surface area obtained from PlanetScope, Basemap, and Planet Fusion. When the PlanetScope surface area date did not correspond exactly to the SkySat dates, we used the closest observation in time, which had a maximum difference of three days before or after the SkySat date. In addition, we assessed the uncertainties of PlanetScope, Basemap, and Planet Fusion for different surface area size classes: 0.1-5 ha (n = 50), 5-10 ha (n = 46), and 10-50 ha (n = 48).

Surface Water Area Validation Using SkySat Imagery
The surface area obtained from PlanetScope, Basemap, and Planet Fusion showed great agreement (r 2 ≥ 0.98) with the validation dataset. In addition, PlanetScope had the smallest MAPE (8.09%), followed by Basemap (8.21%) and Planet Fusion (9.17%) ( Figure 3). When splitting the validation surface area observations into different size classes (Table 4), all three image sources presented similar agreement (r 2 ≥ 0.87), and the highest r 2 values were found for surface area observations between 10 and 30 ha (r 2 ≥ 0.95). All three sources had a similar MAPE for observations between 0.1 and 5 ha (~7.55%) and between 10 and 30 ha (~7.98%), while the highest values were found for observations between 5 and 10 ha (~10.27%). fewer clouds), ground-control ratio (defines the image positional accuracy; values closer to 1 mean higher accuracy), and ground-sampling distance in meters.

Surface Water Area Validation Using SkySat Imagery
The surface area obtained from PlanetScope, Basemap, and Planet Fusion showed great agreement (r 2 ≥ 0.98) with the validation dataset. In addition, PlanetScope had the smallest MAPE (8.09%), followed by Basemap (8.21%) and Planet Fusion (9.17%) ( Figure 3). When splitting the validation surface area observations into different size classes (Table 4), all three image sources presented similar agreement (r 2 ≥ 0.87), and the highest r 2 values were found for surface area observations between 10 and 30 ha (r 2 ≥ 0.95). All three sources had a similar MAPE for observations between 0.1 and 5 ha (~7.55%) and between 10 and 30 ha (~7.98%), while the highest values were found for observations between 5 and 10 ha (~10.27%).

Number of Surface Area Observations per Dataset
The number of PlanetScope observations varied throughout the year and varied across different OFRs (Figure 4). The months with the highest number of PlanetScope images were November-December 2020 (~17) and March-April 2021 (~14), while the months with the lowest numbers were July-September 2020 (<10) and February 2021 (<3) ( Figure 4A). In addition, most of the OFRs (~60%) had 80-100 PlanetScope observations per year ( Figure 4B). Basemap images were processed at a daily cadence, and we considered a new Basemap observation every time a new image composite was used. In this regard, the number of Basemap images followed a similar pattern found for PlanetScope; however, the mean number of Basemap observations per month was higher than of PlanetScope in 10 out of the 12 months analyzed, and most of the OFRs (~75%) had 90-120 Basemap observations per year ( Figure 4A,B).

Planet Fusion Comparison with PlanetScope
We found a high agreement (r 2 ≥ 0.90) for the same-day surface area pairwise comparisons between Planet Fusion and PlanetScope for all size classes ( Figure 5). MAPE decreased as observations increased in size, and the highest MAPE values were found for synthetic, mixed, and real for all size classes ( Figure 5). The number of pairwise comparisons for real was higher than mixed and synthetic, to a large extent (~60%). This finding is somewhat expected, as the Planet Fusion algorithm uses PlanetScope images as an input to generate daily Planet Fusion imagery. Planet Fusion images were derived from real and synthetic pixel values, and the number of real and synthetic observations varied throughout the year ( Figure 4C). The number of images derived from real pixels reached its peak (~13-15) between November and December 2020, and the lowest numbers were found in February 2021 (~2) and between May and June 2021 (<5). In general, most of the OFRs had~80-100 real observations per year. The number of mixed images (i.e., composed by real and synthetic pixels) tended to be < 10 for all months, and most of the OFRs had <50 mixed observations per year ( Figure 4C,D). The number of synthetic images is higher than real and mixed observations for all months of the year, and the highest values (~22-26) occurred in July 2021 and May-June 2020, with the lowest values between November and December 2020 (~14-16) ( Figure 4C). In addition, most of the OFRs had~250-260 synthetic observations per year ( Figure 4D).

Planet Fusion Comparison with PlanetScope
We found a high agreement (r 2 ≥ 0.90) for the same-day surface area pairwise comparisons between Planet Fusion and PlanetScope for all size classes ( Figure 5). MAPE decreased as observations increased in size, and the highest MAPE values were found for synthetic, mixed, and real for all size classes ( Figure 5). The number of pairwise comparisons for real was higher than mixed and synthetic, to a large extent (~60%). This finding is somewhat expected, as the Planet Fusion algorithm uses PlanetScope images as an input to generate daily Planet Fusion imagery.

Monthly Comparisons between Basemap and Planet Fusion with PlanetScope
When comparing each OFR surface area time series derived from Basemap and Planet Fusion with PlanetScope, for both datasets, most of the OFRs (63% and 61% for

Monthly Comparisons between Basemap and Planet Fusion with PlanetScope
When comparing each OFR surface area time series derived from Basemap and Planet Fusion with PlanetScope, for both datasets, most of the OFRs (63% and 61% for Basemap and Planet Fusion, respectively) showed good agreement with r 2 ≥ 0.55, and 74% and 70% of the OFRs presented small uncertainties with MAPE <5% (Figure 6). The mean monthly percent error-calculated by comparing Basemap and Planet Fusion with PlanetScope-for Basemap and Planet Fusion varied between − 2.45-1.48% and between − 3.36%-1.66% for 0.1-5 ha, between − 2.88-1.11% and between −3.56-0.51% for 5-10 ha, and between − 2.23-0.53% and − 3.13-0.76% for 10-30 ha. These values were stable throughout the year (Figure 7). The percent error variability decreased as the surface area observations increased in size, and the observations between 10 and 30 ha had the least variability. In addition, Planet Fusion presented smaller percent error variability when compared to Basemap for all size classes (Figure 7). The highest MAPE values for Basemap (4.73%) and Planet Fusion (5.80%) were found for observations between 0.1 and 5 ha, and the MAPE was <4.40% for all months for both Basemap and Planet Fusion for observations between 5 and 10 ha and 10 and 30 ha, respectively. This indicates that even when there are fewer PlanetScope images available to generate Basemap and Planet Fusion due to clouds, shadow, and haze, both products tend to have surface area uncertainties <5%. The mean monthly percent error-calculated by comparing Basemap and Planet Fusion with PlanetScope-for Basemap and Planet Fusion varied between −2.45-1.48% and between −3.36%-1.66% for 0.1-5 ha, between −2.88-1.11% and between −3.56-0.51% for 5-10 ha, and between −2.23-0.53% and −3.13-0.76% for 10-30 ha. These values were stable throughout the year (Figure 7). The percent error variability decreased as the surface area observations increased in size, and the observations between 10 and 30 ha had the least variability. In addition, Planet Fusion presented smaller percent error variability when compared to Basemap for all size classes (Figure 7). The highest MAPE values for Basemap (4.73%) and Planet Fusion (5.80%) were found for observations between 0.1 and 5 ha, and the MAPE was <4.40% for all months for both Basemap and Planet Fusion for observations between 5 and 10 ha and 10 and 30 ha, respectively. This indicates that even when there are fewer PlanetScope images available to generate Basemap and Planet Fusion due to clouds, shadow, and haze, both products tend to have surface area uncertainties <5%. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 21

OFR Surface area Time Series
We selected six OFRs (Table 2) to illustrate the surface area time series derived from PlanetScope, Basemap, and Planet Fusion (Figure 8). The surface area time series show that different OFRs have different surface area change patterns. In general, the OFR surface area decreased between 07/20 and 11/20 (e.g., Figure 8, OFRs A-D), period of the year when farmers are irrigating their crops [23], and it increased between 01/21 and 05/21, which are the months when the study region receives most of its annual precipitation [23].
When compared to PlanetScope and Basemap, Planet Fusion had a smoother surface area time series with less variability (e.g., Figure 8 OFRs A-D). In addition, the Planet Fusion time series was less affected by the presence of clouds and haze, which can increase or decrease OFR surface area. Even though we used a low cloud-cover threshold (<10%) for PlanetScope, there are several PlanetScope and Basemap images contaminated with cloud shadows and haze (e.g., Figure 8, OFR A, between 08/20 and 09/20), indicating surface area ~20% larger than that of Planet Fusion. Other examples were observed between 07/20 and 08/20 and between 05/21 and 06/21 in Figure 8, OFR B, in which there were no PlanetScope images available and the Basemap shows abrupt changes in surface area-a drop of 20% and 15% for both dates-which were caused by the presence of cloud shadows and haze. In Figure 9, we highlighted the impact of clouds and haze for OFR A (08/16/2020) and OFR B (08/30/2020). For OFR A, PlanetScope and Basemap surface areas were ~20% larger than those of Planet Fusion, which is explained by the misclassification of water on the lower-right corner of the OFR. For OFR B, while the PlanetScope image had a surface area ~13% larger than that of Planet Fusion, the Basemap image indicated a surface area ~14% smaller. These discrepancies are caused by the presence of clouds in the PlanetScope image and haze in the Basemap image.

OFR Surface Area Time Series
We selected six OFRs (Table 2) to illustrate the surface area time series derived from PlanetScope, Basemap, and Planet Fusion (Figure 8). The surface area time series show that different OFRs have different surface area change patterns. In general, the OFR surface area decreased between 20 July and 20 November (e.g., Figure 8, OFRs A-D), period of the year when farmers are irrigating their crops [23], and it increased between 21 January and 21 May, which are the months when the study region receives most of its annual precipitation [23].
When compared to PlanetScope and Basemap, Planet Fusion had a smoother surface area time series with less variability (e.g., Figure 8 OFRs A-D). In addition, the Planet Fusion time series was less affected by the presence of clouds and haze, which can increase or decrease OFR surface area. Even though we used a low cloud-cover threshold (<10%) for PlanetScope, there are several PlanetScope and Basemap images contaminated with cloud shadows and haze (e.g., Figure 8, OFR A, between 20 August and 20 September), indicating surface area~20% larger than that of Planet Fusion. Other examples were observed between 20 July and 20 August and between 21 May and 21 June in Figure 8, OFR B, in which there were no PlanetScope images available and the Basemap shows abrupt changes in surface area-a drop of 20% and 15% for both dates-which were caused by the presence of cloud shadows and haze. In Figure 9, we highlighted the impact of clouds and haze for OFR A (16 August 2020) and OFR B (30 August 2020). For OFR A, PlanetScope and Basemap surface areas were~20% larger than those of Planet Fusion, which is explained by the misclassification of water on the lower-right corner of the OFR. For OFR B, while the PlanetScope image had a surface area~13% larger than that of Planet Fusion, the Basemap image indicated a surface area~14% smaller. These discrepancies are caused by the presence of clouds in the PlanetScope image and haze in the Basemap image.    Table 2) surface-area time series derived from PlanetScope, Basemap, and Planet Fusion and OFR shapefiles overlaid on high-resolution Google Satellite imagery.  Table 2) surface-area time series derived from PlanetScope, Basemap, and Planet Fusion and OFR shapefiles overlaid on high-resolution Google Satellite imagery.

Figure 9.
OFRs A and B (see Table 2) PlanetScope, Basemap, and Planet Fusion false-color composites (blue: red, green: green, and red: NIR) and the surface-water classification for 16 August 2020 (OFR A) and 30 August 2020 (OFR B).
OFR surface water classification is impacted by the OFRs' environmental conditions and their shape geometry. OFRs with complex geometries (e.g., not circular or square and shapes with a large number of edges) tend to have higher surface area classification uncertainties [33]. For example, Figure 8, OFR D, shows a multi-part OFR that may not have all parts inundated at the same time, which can explain part of the variability in the surface area time series for PlanetScope, Basemap, and Planet Fusion. The surface area time series from OFR E and OFR F (Figure 8) are influenced by the presence of vegetation within the OFRs. The presence of vegetation impacts surface water classification [5,33,55], leading to noisy surface area time series and abrupt changes (e.g., OFR E between 20 September and 21 January). In addition, the high variability in surface area for OFRs E and F is related to the presence of adjacent water bodies, which can inundate during flood events and contribute to changes on OFR boundary limits. We highlighted the impact of vegetation on the OFR E time series for two different occasions: 14 July 2020 and 16 October 2020 ( Figure 10). During the first occasion, PlanetScope and Basemap indicated surface area (~9.5 ha) 95% greater than Planet Fusion (0.5 ha); on the second occasion, a contrasting scenario in which Planet Fusion surface area (12.25 ha) was 86% higher than that of PlanetScope and Basemap (~2 ha). These results shed light on the importance of assessing the OFR environmental conditions and how they impact the OFR surface area time series before employing these datasets to monitor surface area changes.  Table 2) PlanetScope, Basemap, and Planet Fusion false-color composites (blue: red, green: green, and red: NIR) and the surface-water classification for 08/16/2020 (OFR A) and 08/30/2020 (OFR B).  Table 2) PlanetScope, Basemap, and Planet Fusion false-color composites (blue: red, green: green, and red: NIR) and the surface water classification for 07/14/2020 and 10/16/2020. Figure 10. OFR E (see Table 2) PlanetScope, Basemap, and Planet Fusion false-color composites (blue: red, green: green, and red: NIR) and the surface water classification for 14 July 2020 and 16 October 2020.

Discussion
The surface area validation carried out using multiple SkySat imagery showed that the methodology used to classify OFR surface area performed well for PlanetScope, Basemap, and Planet Fusion, with high agreement r 2 ≥ 0.87 and MAPE between 7.05% and 10.08% for all image sources and all size classes (Table 4). Comparisons between Basemap and Planet Fusion with PlanetScope highlighted that most of the OFRs had good agreement with 61% of the OFRs with r 2 ≥ 0.55, and small uncertainties with 70% of the OFRs with MAPE < 5% ( Figure 6). Basemap and Planet Fusion presented similar monthly mean percent error (~−3-3%) and MAPE (~2.20-5.80%) throughout the year (Figure 7). In addition, percent error variability and MAPE decreased for the larger surface area observations (Figure 7). The highest monthly MAPE (5.80%) was found for Planet Fusion for observations between 0.1 and 5 ha, and the MAPE was ≤4.40% for Basemap and Planet Fusion for observations between 5 and 10 ha and between 10 and 30 ha. Furthermore, when analyzing the three Planet Fusion data categories (i.e., real, mixed, and synthetic), the greatest uncertainties were found for the synthetic images (MAPE~5%), followed by mixed (MAPE~4%) and real (MAPE~3%) ( Figure 5). These findings indicate that Basemap and Planet Fusion images can be employed to monitor OFRs with uncertainties < 10% when the sources are compared to the validation dataset and with uncertainties < 5% when compared to PlanetScope. However, time series obtained from Basemap and Planet Fusion can be highly variable ( Figure 8E,F), as surface water classification can be impacted by the size of water bodies (Table 4, Figures 2, 4 and 6) and the environment in which OFRs are located (e.g., presence of vegetation within the OFRs; Figure 8, OFRs D-F).
The number of cloud-free observations offered by Basemap and Planet Fusion enlightens the potential of these datasets to monitor OFR surface area changes (Figure 4). Both datasets pose advantages when compared to a single sensor approach-employing Plan-etScope alone (Figure 4), or other sensors, for example, Landsat [23,31,56], Sentinel 1 [57], and Sentinel-2 [58-60]-or a multi-sensor approach [4,33,34]. Briefly, the use of a single sensor is limited to a few observations a month, and in some periods of the year in eastern Arkansas, there could be weeks without a cloud-free image [33]. Although the number of observations is improved when employing a multi-sensor approach, daily to sub-weekly monitoring is not attainable unless an assimilation algorithm [33] is implemented. In addition, when implementing a multi-sensor approach, it is necessary to acquire the data from different platforms (e.g., Planet Explore, Sentinel Hub, and Google Earth Engine), which can be time-consuming and a limiting factor if it is necessary to process, download, and move the satellite imagery across multiple platforms. In this study, we demonstrated that Basemap and Planet Fusion imagery processing can be done entirely in the cloud environment by leveraging the integration of Planet's Platform, Google Cloud Storage, and Google Earth Engine. This integration allows for swift analysis, and it can be used for other study regions without the need to acquire data from multiple platforms.
Daily OFR surface area time series derived from Basemap and Planet Fusion revealed important differences between the two datasets. In general, Basemap had higher surface area variability, and it was more susceptible to the presence of cloud shadows and haze when compared to Planet Fusion, which had a smoother time series with less variability and fewer abrupt changes throughout the year (Figure 8). The Planet Fusion algorithm combines data from multiple satellites to establish a baseline of OFR surface area time series by filling gaps with synthetic pixels. Nonetheless, the smoothing effect should be interpreted cautiously, as some changes in the time series due to large rainfall events or frequent irrigation activities may be smoothed out. This is especially relevant for the periods of the year when there are more synthetic observations ( Figure 4) and the uncertainties in surface area tend to be higher ( Figure 5). Additionally, because Planet Fusion is based on a robust algorithm that uses data from various satellites, this dataset requires more image processing steps and higher computing power when compared to Basemap, which is generated at a faster speed with lower processing costs. Meanwhile, the Basemap time series may contain a "stair-step" effect caused by repeated observations when the Basemap scene composition was kept constant due to a lack of new cloud-free scenes (e.g., Figure 8, OFR D, early March 2021). By keeping the same image composition, the Basemap algorithm avoids generating synthetic pixel values while still providing a cloud-free observation. Nonetheless, it is important to keep in mind that there could be scenarios (e.g., when there is a lack of a new cloud-free scene for weeks or more) in which the Basemap may have the same number of observations as PlanetScope, hence decreasing its monitoring capabilities.
Our findings have important implications to future hydrological studies that aim to monitor small water bodies at large scale and high temporal frequency. For the OFRs in eastern Arkansas, the Basemap and Planet Fusion surface area time series helped unravel sub-weekly changes in OFR surface area, as well as yearly seasonality (Figure 8). OFRs surface area changes are pivotal information for calculation of OFR water volume inflows and outflows. This can be achieved by combining the surface area time series with the area-volume equations (e.g., hypsometry), which are derived using the OFRs' geometric shape and depth [56,61,62]. Estimating OFRs volume change helps bridge one of the key limitations when modeling the cumulative impacts of OFRs on surface hydrology, as OFR water volume change is commonly assumed to be equal to all OFRs located in a watershed [24,25,63]. In addition, as the number of OFRs is projected to increase globally [24,27], understanding the impact of OFRs on surface hydrology is pivotal when seeking indicators to determine the optimal spatial distribution and number of OFRs, as well as their storage capacities and water management plans aiming to mitigate downstream impacts. Beyond implications to hydrological studies, we demonstrated that Basemap and Planet Fusion can be used to monitor surface area changes for a network of OFRs (Figure 8). This information can be used by regulatory agencies to create water status reports to improve regional water management and water use efficiency. These reports would be especially relevant during the dry critical period of the year when farmers are frequently irrigating.

Known Issues and Limitations
We applied Basemap and Planet Fusion imagery for a one-year analysis. More research is necessary to assess the performance of these datasets for a longer study period (e.g., including periods of prolonged droughts~3-5 years) and in other study regions-for example, in Southern India, where OFRs are common [4], and where there is a monsoon climate in which there could be weeks without a clear-sky image [64]. In addition, the validation of this study was conducted using cloud-free SkySat imagery; therefore, there is still a need to further evaluate the performance of both datasets under cloudy conditions. However, this will require extensive field work, including visiting multiple OFRs on cloudy days, which imposes several challenges, as most of the OFRs in eastern Arkansas are located on private properties. Furthermore, we assumed that OFR surface area would vary within known and limited boundaries (i.e., OFR shapefile buffered to 20 m). However, different results might be obtained if the Basemap and Planet Fusion images are used to monitor water bodies that frequently change their boundaries-water impoundments that are located close to water streams and rivers that flood frequently, impacting the edges of water bodies. Lastly, although we calculated the uncertainties introduced by Basemap and Planet Fusion, when using these datasets for monitoring purposes, it would be helpful to have an estimated uncertainty accompanying every surface area observation. For instance, whenever there are repeated observations by the Basemap or continuous synthetic observations from Planet Fusion, the uncertainties from these images will be higher; however, as of now, we cannot estimate an observation based uncertainty.

Conclusions
We presented a novel application of Basemap and Planet Fusion analysis ready datasets to monitor sub-weekly OFRs surface area changes. We tested both datasets to monitor 340 OFRs of different sizes, and we found that these datasets can be employed to monitor OFRs with uncertainties < 10% when compared to an independent validation dataset and with uncertainties < 5% when compared to PlanetScope imagery. While Basemap had higher surface area variability and it was more susceptible to the presence of cloud shadows and haze, Planet Fusion had a smoother time series with less variability and fewer abrupt changes throughout the year. Given that the surface area classification can be impacted by the OFR environmental conditions (e.g., presence of vegetation inside the OFR), therefore limiting the use of these datasets, we recommend assessing the OFRs' surface area time series before employing them for monitoring purposes. As the number of OFRs is expected to increase globally, the use of these datasets is of great importance to understanding OFR sub-weekly, seasonal and inter-annual surface area changes, and to improving freshwater management by allowing better assessment and management of OFRs.