Multiscale Radiometric Stability Analysis of Water Bodies in Multispectral Remote Sensing Imagery

TL;DR

This paper shows that using flux-conserving resampling helps maintain accurate radiometric data when combining multispectral remote sensing images of water bodies at different scales.

Contribution

The study introduces a radiometric stability-based framework to evaluate resampling methods for multiscale aquatic remote sensing.

Findings

Radiometric consistency decreases with increasing spatial scale in aquatic remote sensing.

Flux-conserving resampling maintains higher radiometric stability and preserves water radiance characteristics.

The proposed framework supports reliable multi-source data fusion and quantitative inversion in aquatic remote sensing.

Abstract

What are the main findings?
Spatial resampling prior to data fusion can introduce substantial radiometric distortions in aquatic remote sensing, particularly at coarse spatial resolutions.Flux-conserving resampling consistently preserves radiometric stability across scales.

Spatial resampling prior to data fusion can introduce substantial radiometric distortions in aquatic remote sensing, particularly at coarse spatial resolutions.

Flux-conserving resampling consistently preserves radiometric stability across scales.

What are the implications of the main findings?
Enforcing flux conservation during resampling is critical for maintaining the integrity of weak aquatic signals and improving the reliability of water quality retrievals.The results provide a physically grounded reference for multi-source aquatic data fusion and cross-calibration of sensors with differing spatial resolutions.

Enforcing flux conservation during resampling is critical for maintaining the integrity of weak aquatic signals and improving the reliability of water quality retrievals.

The results provide a physically grounded reference for multi-source aquatic data fusion and cross-calibration of sensors with differing spatial resolutions.

In remote sensing, multi-sensor data fusion enhances environmental monitoring by integrating complementary observations. A critical step in this integration is spatial resampling to a common scale. Although often regarded as a routine preprocessing operation, resampling can become a significant source of radiometric uncertainty, systematically altering scene radiance during scale transformation, especially in heterogeneous aquatic environments. In this study, we evaluate resampling-induced radiometric uncertainty and assess the physical advantages of flux-conserving resampling in multi-scale aquatic remote sensing. Using the radiometrically stable Landsat 8 OLI sensor as a reference platform, this study develops a radiometric stability–based framework to evaluate multi-scale resampling methods. Radiometric consistency in the visible bands was first evaluated using a Rayleigh scattering calibration, allowing a systematic comparison of four resampling methods across multiple spatial scales. Normalized water-leaving radiance was then retrieved using the Satellite Signal in the Solar Spectrum (6S) radiative transfer model and validated against in situ AERONET-OC measurements. Our results indicate that radiometric consistency decreases with increasing scale, while flux-conserving resampling maintains higher stability and preserves the spatiotemporal characteristics of water radiance. These findings highlight the importance of flux-conserving resampling for multi-scale radiometric fidelity and establish the proposed framework as a reference for reliable multi-source data fusion and quantitative inversion in aquatic remote sensing and beyond.

Figures

Full text

1. Introduction

Satellite remote sensing plays a crucial role in monitoring Earth’s systems [1,2]. Among its applications, multispectral imagery serves as a fundamental tool for quantitative environmental analysis [3]. This capability is particularly evident in ocean color remote sensing, where spectral reflectance supports large-scale and cost-effective monitoring of water quality [4]. However, dedicated ocean color sensors often suffer from coarse spatial resolution, which limits their ability to monitor dynamic processes in coastal and inland waters [5]. Satellites such as Landsat provide relatively higher spatial resolution (30 m), stable calibration, and global coverage. These features make them valuable data sources for aquatic studies [6,7,8]. Nevertheless, passive optical sensors face inherent trade-offs among spatial, spectral, and temporal resolution because photon collection limits and signal-to-noise constraints [9]. Consequently, observations from a single satellite platform are often insufficient for comprehensively characterizing aquatic optical properties, particularly in optically heterogeneous environments and across varying spatial scales.

To overcome these limitations, multi-source data fusion has become an increasingly common strategy to integrate complementary sensors and achieve enhanced spatio-spectral resolution [10]. Extensive research has addressed multi-source data fusion, particularly in relation to spatiotemporal integration [11] and the harmonization of heterogeneous datasets with differing spatial resolutions [12,13]. Among the technical considerations for multi-source fusion, a fundamental requirement is spatial scale conversion, which involves resampling imagery from different sensors onto a common grid [14]. Despite its ubiquitous use, resampling is often treated as a purely technical or geometric operation. In practice, however, it can substantially alter spectral information and compromise pixel-level radiometric integrity, particularly in optically heterogeneous aquatic environments [15]. Previous multi-scale analyses of water bodies have shown that the choice of resampling algorithm affects both the preservation of spatial detail and the consistency of spectral and radiometric signals [16]. Similar effects have also been reported in high-resolution UAV studies, where spatial resolution and resampling strategies were shown to influence classification accuracy and object discrimination [17]. These effects can directly propagate into downstream quantitative products, such as water-leaving radiance, yet they remain insufficiently examined in aquatic remote sensing workflows.

In practice, most studies rely on computationally efficient interpolation-based methods, such as bilinear interpolation, cubic convolution, or Lanczos resampling [18,19]. Although widely adopted, these conventional interpolation methods do not preserve radiative flux during spatial scale transformation and may distort scene-integrated radiance, leading to systematic radiometric and spectral errors. For instance, Tian et al. showed that as spatial resolution decreases, traditional interpolation tends to reduce the depth and area of key spectral features, producing quantifiable spectral distortions compared with more physically informed resampling methods [20]. Importantly, such errors are not limited to traditional image processing workflows: they can also propagate into downstream tasks, such as the quality of synthesized data used in deep learning–based super-resolution training, thereby affecting analytical performance [21]. In contrast, flux-conserving resampling is grounded in the physical principle of energy conservation and preserves total radiance during spatial scale transformation, mitigating these errors. Its effectiveness has been demonstrated in marine remote sensing for harmonizing observations from heterogeneous sensors within virtual constellations [22], as well as in atmospheric remote sensing, where flux-conserving downscaling improves trace gas retrieval accuracy by maintaining radiometric consistency [23].

Despite the promising potential of flux-conserving methods across various remote sensing domains, a systematic evaluation of resampling methods specifically designed for multi-scale aquatic remote sensing remains understudied. In particular, there has been limited effort to evaluate how different resampling methods affect radiometric stability across spatial scales or how such effects propagate into physically meaningful end products, such as normalized water-leaving radiance.

In this study, we focus on the methodological role of spatial resampling as a non-negligible source of radiometric variability in water color remote sensing. Using Landsat-8 OLI data, this study assesses four commonly applied resampling methods, namely bilinear interpolation, cubic convolution, Lanczos interpolation, and flux-conserving resampling, across multiple spatial scales (500 m, 1000 m, and 2000 m). Radiometric stability in the visible bands is first examined using a calibration framework based on Rayleigh scattering. The influence of resampling on normalized water-leaving radiance is then quantified with the 6S radiative transfer model and validated against in situ AERONET-OC observations.

This study demonstrates that spatial resampling is a critical methodological component rather than a simple technical preprocessing step in multi-scale aquatic remote sensing, by directly linking the choice of resampling strategy to radiometric consistency and water color retrieval accuracy. The results highlight the advantages of flux-conserving resampling for maintaining radiometric fidelity across scales and provide a practical framework for implementing reliable, radiometrically consistent multi-source data fusion in aquatic environments.

2. Study Area and Data

2.1. Overview of the Study Area

In this study, the nearshore waters of the eastern Gulf of Thailand (7.5–9.5° N, 100.0–102.0° E; Figure 1) were selected as the primary study area. This region, situated between the Indochina and Malay Peninsulas, comprises a typical tropical semi-enclosed shallow bay with an average depth of 40–60 m and a gentle bathymetric slope, characteristic of continental shelf seas [24]. The strong influence of the Chao Phraya River discharge generates a pronounced optical gradient from highly turbid nearshore waters to clearer offshore waters.

The region is influenced by the Asian monsoon, which drives seasonal variations in hydrographic and atmospheric conditions. Under the northeast monsoon during the dry season (November–February), precipitation and river discharge are reduced, and atmospheric conditions are relatively stable with low aerosol loading, resulting in more homogeneous water optical properties. These relatively stable atmospheric and hydrographic conditions reduce uncertainties in atmospheric correction, thereby providing an ideal natural setting for controlled validation of radiometric calibration and water color remote sensing algorithms. In contrast, the wet season is characterized by higher rainfall and stronger river inflow, which increase water turbidity and variability in atmospheric aerosol conditions. Therefore, the study focused on dry-season conditions to ensure controlled and reliable radiometric validation, thereby minimizing the confounding effects of atmospheric complexity.

2.2. Data Sources

The Landsat-8 Operational Land Imager (OLI) scene used in this study was acquired on 15 December 2015 (Path/Row: 128/54). Landsat-8 operates in a sun-synchronous orbit at an altitude of 705 km with a 16-day revisit cycle [25]. The OLI sensor provides multispectral data at 30 m spatial resolution (15 m for the panchromatic band) across a 185 km swath, featuring improved radiometric calibration accuracy and the signal-to-noise ratio compared to previous Landsat missions [26]. Although the scene was acquired in 2015, it was selected because it captured optimal atmospheric and water conditions during the dry season, providing a reliable reference for evaluating the resampling algorithms. The spectral characteristics of the OLI bands are summarized in Table 1, including the coastal aerosol band (443 nm), which is particularly important for aquatic applications.

To evaluate the performance of the resampling methods across different water conditions, three additional AERONET-OC sites were included: MVCO (41.3° N, 70.6° W), Helsinki Lighthouse (59.9° N, 24.9° E), and Lucinda (18.5° S, 146.4° E). Landsat-8 scenes for the auxiliary sites were selected primarily based on cloud cover (<10%) without seasonal restrictions. All scenes, including those from the primary study area in the eastern Gulf of Thailand, were used for water-leaving radiance retrieval to assess the performance of the resampling methods. Scene selection followed strict criteria to ensure data quality, and detailed information for all matched scenes was provided in Table 2.

The primary validation was conducted using the Landsat-8 scene matched to the GOT-Seaprism site (9.3° N, 101.4° E) of the Aerosol Robotic Network for Ocean Color (AERONET-OC) [27], which served as the main reference for radiometric stability analysis. Satellite observations were temporally matched with in situ measurements by averaging all available data acquired within a ±30 min window centered on the satellite overpass time. Data from the three auxiliary sites were used to evaluate the performance of the resampling algorithms across diverse aquatic environments. The geographic distribution of the primary study area and all validation sites is shown in Figure 2.

3. Methods

3.1. Radiometric Calibration and Atmospheric Correction Models

3.1.1. Radiometric Calibration Model

In this study, Level-1 Landsat-8 OLI data were used as a reference for evaluating radiometric errors introduced by spatial resampling. The OLI instrument has a rigorous radiometric calibration framework, including pre-launch absolute calibration using an integrating sphere and a solar simulator traceable to the National Institute of Standards and Technology (NIST), as well as continuous on-orbit monitoring via a dual solar diffuser system [28]. Such radiometric stability is critical for this study, as it minimizes sensor-induced uncertainties and allows observed radiometric differences to be primarily attributed to spatial resampling effects rather than instrumental artifacts. Top of atmosphere (TOA) radiance was computed from digital number (DN) values using the standard calibration coefficients:

$$[eqn]$$

where $[eqn]$ denotes the TOA radiance, DN represents the digital number of a pixel, and gain and offset are the radiometric calibration coefficients corresponding to the sensor’s gain and bias. The TOA radiance was subsequently converted to TOA reflectance using the following expression:

$$[eqn]$$

where $[eqn]$ converts the directional radiance $[eqn]$ to hemispherical reflectance, $[eqn]$ is the Earth-Sun distance correction factor for orbital eccentricity, $[eqn]$ is the mean exoatmospheric solar irradiance, and $[eqn]$ corrects for the projection of irradiance onto a horizontal surface based on the solar zenith angle $[eqn]$ .

3.1.2. Atmospheric Correction Model

A controlled experiment to isolate spatial resampling effects was conducted within the 6S radiative transfer framework using its deterministic ocean surface model under fixed atmospheric conditions. The ocean surface reflectance $[eqn]$ was calculated using the standard 6S formulation [29]:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ are the solar and sensor azimuth angles, and $[eqn]$ is the relative azimuth angle. $[eqn]$ is whitecap reflectance, $[eqn]$ is sun glint, $[eqn]$ is water-leaving reflectance, and W is whitecap coverage. The ocean surface reflectance was treated as a fixed input, while top of atmosphere reflectance was simulated from this surface model using the 6S radiative transfer calculations, ensuring that variations across spatial resolutions were attributable to pixel aggregation.

3.1.3. Lookup Table Construction

The lookup table (LUT) was constructed using the 6S radiative transfer model to support calibration validation. By precomputing reflectance values across a range of atmospheric and geometric conditions, the LUT holds these factors constant, ensuring that observed radiometric differences between images at different resolutions primarily reflect the effects of spatial resampling.

For the GOT-Seaprism site, the LUT simulated water-leaving reflectance across the full range of expected conditions, including variations in solar and sensor zenith angles, relative azimuth angles, aerosol loading, and sea surface roughness (represented by wind speed). The LUT was populated with realistic data, including Aerosol Optical Depth (AOD) from MODIS and wind speed, water vapor, and ozone from the ERA5 reanalysis. The ranges of key variables such as AOD and wind speed were determined from statistical analysis over the study area. The chosen discretization steps were set to balance computational efficiency and interpolation accuracy, with step sizes sufficiently fine to capture typical atmospheric and geometric variability. This choice ensures that interpolation-induced errors remain limited, as further quantified in Section 5.1.2. Water vapor (3.75 g $[eqn]$ cm^−2^) and ozone (547.46 DU) were set to representative, stable values. In this study, these two atmospheric constituents were held constant to ensure consistency across resampling experiments and to isolate the effect of spatial processing. Although both water vapor and ozone may exhibit spatial variability at larger spatial scales, their horizontal gradients within the selected marine study area are relatively limited under the atmospheric conditions considered. Treating them as spatially uniform, therefore, enables a controlled comparison of radiometric differences induced by resampling, without introducing additional variability from atmospheric heterogeneity. Reflectance for a specific scene was then obtained via linear interpolation within the fixed table. The distributions of LUT input parameters are summarized in Figure 3 and Table 3.

3.2. Multiscale Resampling Algorithm

Spatial resampling is a critical step in spatial scale transformation in remote sensing imagery. Ocean color signals are particularly vulnerable to distortions during this process due to their weak water-leaving radiance and high sensitivity to atmospheric and sub-pixel effects. Conventional interpolation methods, which do not conserve radiometric flux, can further exacerbate these distortions, leading to signal attenuation.

In this study, the resampling method was treated as the only controlled variable to quantitatively evaluate its impact on radiometric consistency and ocean color retrieval accuracy. To isolate resampling-induced effects, observation geometry, including BRDF and solar–viewing parameters, atmospheric correction procedures, and masking strategies were kept identical across all experiments, while only the pixel resampling algorithm was varied. Under this controlled framework, a flux-conserving (FC) resampling method was applied, while bilinear, cubic convolution, and Lanczos interpolation were used for comparison in multi-scale experiments.

3.2.1. Flux-Conserving Resampling

Flux-conserving resampling is an image resampling method based on the Sutherland–Hodgman polygon clipping algorithm [30]. It is grounded in the flux-conserving scheme originally proposed by McGlynn [31]. Unlike conventional resampling techniques, it incorporates the physical principle of total flux conservation, ensuring that the total radiometric energy within the observation area remains unchanged during pixel-scale transformation.

In remote sensing, total flux refers to the sum of the radiance values of each pixel multiplied by its corresponding spatial area, representing the total radiometric energy contribution of the entire region. The flux-conserving approach calculates the overlap area between source and target pixels and applies area-weighted accumulation of the original radiance values to the target grid cells. This process preserves radiometric flux and ensures physical consistency during scale conversion.

As illustrated in Figure 4, the source pixel grid serves as the clipping window. Each resampled target pixel $[eqn]$ is assumed to overlap spatially with its nine surrounding source pixels ( $[eqn]$ ). Each overlapping region is decomposed into triangles using a concave polygon clipping algorithm, and the total area of the overlap is calculated. The overlap fraction $[eqn]$ for each source pixel $[eqn]$ relative to $[eqn]$ is then determined. Based on these overlap coefficients, the flux contribution $[eqn]$ of each source pixel to $[eqn]$ is computed, yielding the resampled pixel value $[eqn]$ as the weighted sum of all contributions:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The flux-conserving algorithm was implemented in Python 3.9 for this study, adapting the core area-weighting scheme. Our implementation supports multi-band imagery and applies proportional clipping to edge pixels, such that only the valid overlapping portions of input pixels contribute to output pixels at boundaries, ensuring radiometric consistency across spatial scale transformations.

3.2.2. Alternative Resampling Methods

Bilinear interpolation

In bilinear interpolation, the value of a target pixel is estimated from the four nearest source pixels (A, B, C, and D) through distance-weighted averaging. As shown in Figure 5a, linear interpolation is first carried out along the x-axis between adjacent pixels, followed by interpolation along the y-axis, yielding the final value at point P (x, y).

2.Cubic convolution interpolation

Cubic convolution interpolation calculates the value at point P (i + u, j + v) using the 16 pixels within a 4 × 4 neighborhood. Each contribution is weighted according to a cubic convolution function, as illustrated in Figure 5b.

3.Lanczos interpolation

Lanczos interpolation estimates the target pixel value by weighted averaging within a 2a × 2a neighborhood, where a = 4 in this study. As shown in Figure 5c, the weighting is determined by the Lanczos kernel function:

$$[eqn]$$

The kernel function features a central lobe of width a, flanked by $[eqn]$ alternating negative and positive side lobes. $[eqn]$ denotes the original pixel values. The two-dimensional interpolation for a target pixel is computed as the product of two one-dimensional Lanczos kernels:

$$[eqn]$$

3.3. Water Pixel Selection Method

To evaluate the impact of spatial resampling on radiometric stability, we developed a workflow for selecting clean water pixels. Throughout the process, steps such as land masking, cloud detection, and glint filtering ensured stable atmospheric and surface water conditions, so that variations in water pixel radiance could be attributed primarily to the resampling method applied. Finally, water regions were extracted using a dual-band thresholding approach based on the green and near-infrared (NIR) bands.

3.3.1. Land and Cloud Mask

Land and cloud pixels were removed prior to water pixel selection to ensure radiometrically homogeneous conditions. Owing to their distinct spectral characteristics, water pixels typically exhibit low values of the normalized difference vegetation index (NDVI), in contrast to the higher values associated with vegetation and bare soil. Therefore, a land mask was created by applying a threshold of NDVI > 0.08, which was empirically determined through iterative testing to ensure stable land–water separation under local conditions:

$$[eqn]$$

Here, NIR and RED represent the reflectance in the near-infrared and red bands. Optimal NDVI thresholds may vary in other regions and conditions and should be adjusted accordingly [32].

Cloud pixels are characterized by high reflectance in the SWIR1 band (B6) [33]. Based on analysis of the imagery characteristics in the study area and multiple tests, an empirical threshold was adopted for initial cloud detection:

$$[eqn]$$

Under the conditions considered in this study, this threshold reliably identifies cloud pixels while minimizing misclassification of water areas. A morphological closing operation (dilation followed by erosion) was subsequently applied to refine cloud boundaries and generate a spatially coherent cloud mask.

3.3.2. Sun-Glint Filtering

To mitigate specular reflection effects, pixels affected by sun glint were excluded using a geometry-based criterion. The glint angle was computed from solar and sensor zenith and azimuth angles provided in the Landsat-8 metadata. Pixels with glint angles ≤ 40° were removed, following established thresholds reported in previous studies [34].

$$[eqn]$$

where $[eqn]$ is the glint angle, $[eqn]$ and $[eqn]$ are the solar and sensor zenith angles, and $[eqn]$ (in Equation (3)) denotes the relative azimuth angle.

3.3.3. Dual-Band Water Pixel Identification

After land, cloud, and glint removal, clear water pixels were identified using a dual-band thresholding approach applied to TOA reflectance. Pixels were classified as water only when reflectance in both the green and NIR bands satisfied predefined thresholds ( $[eqn]$ < 0.08 and $[eqn]$ < 0.15) [35,36], which ensures suppression of atmospheric path radiance and dark non-water targets.

$$[eqn]$$

Sensitivity tests indicated that moderate variations in threshold values did not affect the relative performance ranking of resampling methods.

3.4. Radiometric Stability Validation Methodology

3.4.1. Rayleigh Scattering Calibration

Calibration Site Selection and Water Pixel Stability Analysis

To provide a spatially homogeneous and spectrally consistent baseline for radiometric validation, candidate open-water calibration regions were identified using a sliding-window approach designed to minimize land, cloud, and adjacency contamination, ensuring that subsequent radiometric comparisons primarily reflect resampling effects. The analysis was conducted in the UTM projection to maintain metric consistency, with a window size of 100 km × 100 km and a 20 km step to provide spatially representative, overlapping samples. Water pixels within each window were identified using the water mask and preprocessing scheme described in Section 3.3.

The selection procedure involved three steps:

• (1)Water coverage screening: to ensure high spatial homogeneity and minimize contamination from land, clouds, or adjacency effects, only windows with water coverage exceeding 90% were retained;
• (2)Spatial connectivity validation: an 8-neighborhood connected component analysis was applied, and windows were kept only if the largest contiguous water body accounted for >90% of the total area;
• (3)Coordinate conversion: centroids of the selected windows were transformed from UTM to WGS84 geographic coordinates for geolocation.

The screening procedure initially yielded 12 candidate windows, from which 4 spectrally homogeneous regions were identified. The window with the highest water coverage among these four was selected as a representative calibration site (Figure 6), with a geographic centroid at 8.95° N, 100.95° E, covering a 100 km × 100 km area and a water coverage of 94.80% at the native 30 m resolution. As shown in Figure 6a, the red box marks the selected calibration site; the black areas indicate regions with no data, and the colored area represents the actual Landsat-8 scene extent.

At the study site, we quantitatively evaluated four spatial resampling methods—bilinear interpolation, cubic convolution, Lanczos, and flux-conserving—for their effects on water pixel spatial consistency at 500, 1000, and 2000 m. Figure 7a reports the proportion of valid water pixels at each resampled resolution relative to the 30 m baseline (0.9480). Figure 7b shows boxplots of water pixel stability; the standard deviation (σ) quantifies variability across scales.

Results show that cubic convolution and Lanczos resampling produce higher proportions of water pixels at coarser resolutions (maximum: 0.9608), primarily due to smoothing of high-frequency textures and boundary blurring. Bilinear interpolation exhibits pronounced instability across scales (σ = 0.0105), indicating limited spatial consistency. In contrast, flux-conserving resampling achieves the highest stability (σ = 0.0026), effectively preserving water pixel structures during scaling. These results demonstrate its superiority for applications sensitive to radiometry.

2.Rayleigh Scattering Radiometric Validation

This approach leverages molecular Rayleigh scattering over the open ocean as an absolute radiometric reference [37]. Validation is performed by comparing Landsat-8 OLI visible-band measurements (483–655 nm) with 6S-simulated values interpolated from precomputed LUTs. The Rayleigh calibration coefficient (A) is defined as follows [38]:

$$[eqn]$$

Here, $[eqn]$ represents the observed TOA reflectance from Landsat-8 OLI, and $[eqn]$ denotes the corresponding 6S-simulated TOA reflectance under identical geometric and atmospheric conditions.

Based on the models and methods described in the preceding sections in Section 3, the complete workflow for calculating this coefficient is illustrated in Figure 8.

3.4.2. Retrieval of Normalized Water-Leaving Radiance

In ocean color remote sensing, the water-leaving radiance ( $[eqn]$ ) is derived by subtracting surface-reflected radiance and atmospheric scattering contributions from the total radiance measured at the top of the atmosphere [39]. This process is expressed as:

$$[eqn]$$

$[eqn]$ is the total radiance at the top of the atmosphere; $[eqn]$ and $[eqn]$ are the diffuse and direct atmospheric transmittances; $[eqn]$ and $[eqn]$ denote the whitecap and glint radiances; $[eqn]$ accounts for Rayleigh scattering and its interactions with aerosols; and $[eqn]$ represents the path radiance due to aerosols.

Atmospheric correction was performed using the 6S radiative transfer model. The procedure involved sequentially removing the contributions of Rayleigh scattering and aerosol path radiance from the top of atmosphere radiance to obtain the sea surface radiance. Subsequently, surface-related contributions, including sun glint and whitecap radiance, were subtracted from the sea surface radiance to obtain the water-leaving radiance at each wavelength.

To facilitate quantitative comparison of aquatic optical properties across different observation scenarios, the water-leaving radiance at each wavelength is converted to normalized water-leaving radiance. The complete retrieval workflow is illustrated in Figure 9, and the normalized water-leaving radiance is computed according to Wang et al. [40] as:

$$[eqn]$$

In this equation, $[eqn]$ denotes the normalized water-leaving radiance; $[eqn]$ is the mean extraterrestrial solar irradiance; $[eqn]$ represents the downwelling irradiance just above the water surface; $[eqn]$ is the Earth–Sun distance correction factor; $[eqn]$ corresponds to the Fresnel reflectance of the sea surface; $[eqn]$ is the downward atmospheric transmittance; and $[eqn]$ accounts for the solar zenith angle.

4. Results

4.1. Multi-Scale Validation of Rayleigh Scattering Radiometric Stability

The visible bands of Landsat-8 OLI (443 nm, 486 nm, 561 nm, and 655 nm) cover the shortwave spectrum spanning blue to red, where Rayleigh scattering is pronounced. These bands are therefore suitable for evaluating the accuracy of Rayleigh scattering models and the stability of radiometric calibration coefficients. Clean water pixels at 500 m, 1000 m, and 2000 m spatial resolutions were extracted as described in Section 3.3. Together with satellite observation geometry and 6S radiative transfer model LUTs, pixel-level top of atmosphere (TOA) apparent reflectance was simulated.

The Rayleigh scattering calibration coefficient A was derived from the ratio of measured TOA reflectance (based on Landsat-8 OLI radiometric calibration) to 6S-simulated calibrated TOA reflectance (Equation (13)). A linear regression model was then established between these two datasets, as shown in Figure 10, Figure 11 and Figure 12. In the figures, the x-axis represents measured TOA reflectance from Landsat-8 OLI, and the y-axis represents the 6S-simulated calibrated TOA reflectance. The scatter points correspond to pixel-level data. The solid line represents the fitted regression line, while the dashed line denotes the 1:1 reference line.

Radiometric calibration performance exhibits marked variations among different resampling methods and spatial resolutions (Figure 10, Figure 11 and Figure 12). At 500 m resolution, all methods closely follow the 1:1 reference line, with flux-conserving (FC) resampling showing the tightest clustering. As pixel size increases to 1000 m, observable pixel-mixing effects reduce calibration accuracy, while coarsening to 2000 m leads to pronounced radiometric smoothing, larger deviations, and increased scatter. Method-dependent differences become more evident at coarser resolutions, particularly in the blue bands, where FC resampling outperforms bilinear interpolation by maintaining superior agreement with reference data. Calibration coefficients (A) display wavelength-dependent variations, with larger changes in shorter bands (B1) and minimal variation in longer bands (B4), reflecting the differential sensitivity of spectral bands to resampling.

These observations are consistent with the underlying principles of each resampling technique. The energy-conserving approach of FC resampling mitigates nonlinear error accumulation during scale transformations, whereas bilinear and cubic convolution are more susceptible to subpixel heterogeneity. Lanczos interpolation, while reducing high-frequency artifacts via its sinc-based kernel, does not fully conserve radiometric energy, leading to moderate performance degradation at coarser resolutions.

Overall, FC resampling consistently outperforms alternative methods, preserving radiometric fidelity and maintaining stable calibration coefficients across all spectral bands, even under coarse-resolution conditions with pronounced mixed-pixel effects.

To quantify these observations, relative errors (δ) were calculated to assess radiometric consistency across spatial resolutions for each spectral band as defined in Equation (16):

$$[eqn]$$

where $[eqn]$ denotes the measured TOA reflectance from Landsat-8 OLI and $[eqn]$ represents the 6S-calibrated TOA reflectance, respectively, consistent with the notation in Equation (13) for the calibration coefficient A.

For each resampling method at the three investigated spatial resolutions, relative errors were computed for all spectral bands. The results are presented in Figure 13 and Table 4. In Figure 13, the boxplots summarize the statistical distribution of errors: the boxes indicate the interquartile range (IQR, 25th–75th percentiles), the whiskers extend to 1.5 × IQR, the horizontal bars represent the minimum and maximum values within this range, and the solid line within each box indicates the median.

As shown in Figure 13 and Table 4, the flux-conserving (FC) resampling method consistently demonstrated superior performance across all investigated spatial resolutions (500–2000 m) and spectral bands (443–655 nm), achieving both the lowest relative errors and the highest stability. For instance, in the 443 nm band, FC maintained an average relative error of 5.71% at 500 m resolution, increasing only moderately to 6.93% at 2000 m. In contrast, conventional interpolation methods exhibited more pronounced error growth across the same resolution range: bilinear (BL) increased from 5.90% to 7.80%, while cubic convolution (CC) rose from 6.74% to 7.17%. Lanczos interpolation showed performance comparable to FC at finer resolutions (e.g., 5.60% at 500 m for 486 nm), but degraded more noticeably at coarser resolutions.

In addition to its overall accuracy, distinct spectral dependencies were observed in resampling accuracy. Shorter wavelength bands (443 nm and 486 nm) consistently exhibited higher mean relative errors (5.26–7.80%) and greater variability (standard deviations of 3.48–5.31%) compared to longer wavelengths (561 nm and 655 nm), likely due to stronger spatial heterogeneity and increased sensitivity to atmospheric scattering. Notably, FC maintained exceptional consistency across the spectrum, with standard deviations remaining below 4.64% even at the coarsest resolution (2000 m). This observed superiority in radiometric stability aligns with the findings of Lin et al. [41] in multi-sensor data fusion, underscoring the generalizability of the flux-conserving principle for remote sensing applications requiring high radiometric fidelity.

4.2. Multi-Scale Water-Leaving Radiance Retrieval Validation

To evaluate the retrieval accuracy and radiometric consistency of the multi-scale methodology, the normalized water-leaving radiance $[eqn]$ was retrieved and validated at four geographically diverse AERONET-OC sites: GOT-Seaprism, MVCO, Helsinki Lighthouse, and Lucinda. Landsat-8 imagery at 30 m resolution was used as the baseline and resampled to coarser spatial resolutions of 500 m, 1000 m, and 2000 m using four different resampling methods. For each scale, clear-water pixels were identified and atmospherically corrected using the 6S radiative transfer model to obtain $[eqn]$ in the visible bands. For validation, the mean retrieved $[eqn]$ values at each scale were calculated to represent the regional radiometric characteristics. These values were then compared with concurrent or nearest-available Level 1.5 measurements from the corresponding AERONET-OC sites, using the temporal matching window defined in Section 2.2. This procedure enables assessment of both scale-dependent retrieval performance and absolute radiometric accuracy.

Figure 14 and Table 5 present a comprehensive evaluation of resampling performance across four geographically distinct aquatic sites, revealing consistent patterns in radiometric behavior. Bilinear interpolation consistently produces higher relative errors, particularly in the red band (655 nm), and these errors increase with coarser spatial resolutions, indicating strong sensitivity to boundary effects in heterogeneous waters. Cubic convolution performs competitively in some cases but shows substantial variability across sites, reflecting its dependence on local radiometric structure. Lanczos resampling achieves moderate accuracy but demonstrates a systematic underestimation of reflectance across most bands.

In contrast, the flux-conserving (FC) method exhibits overall superior radiometric stability, achieving the lowest errors in 32 of the 48 evaluated band–resolution–site combinations. Notably, FC maintains relative errors within ±5% across the blue–green spectral range (443–561 nm) for the majority of sites and resolutions, indicating reliable radiometric preservation across spatial scales. Accuracy degrades slightly at the coarsest resolution (2000 m), as exemplified by the −10.0% error at 561 nm for GOT-Seaprism, reflecting the challenges of extreme pixel aggregation.

In addition, the optimal resampling method occasionally varies with site and spectral band. At the Lucinda site, characterized as clear, optically deep water, Lanczos achieved slightly lower relative errors in the blue band (443 nm) at 500 m resolution. A similar pattern is observed at the Helsinki Lighthouse site (500 m, 443–486 nm), where cubic convolution or Lanczos performs competitively.

Interpolation-based methods also occasionally outperform FC in select red-band cases (655 nm), such as GOT-Seaprism and Lucinda at 500 m resolution, as well as in some coarser-resolution combinations (1000–2000 m) for several sites. These occurrences likely reflect relatively homogeneous radiometric structures with weak spatial gradients, where interpolation-based smoothing can incidentally reduce high-frequency noise, residual atmospheric effects, or edge-induced errors near land–water boundaries.

Despite these minor variations, the differences are small (typically within ±2–3%) and do not alter the overall conclusion that FC provides robust radiometric stability across diverse aquatic regimes.

Finally, the persistent negative bias in the red band (655 nm) observed across multiple sites is likely due to residual errors in atmospheric correction rather than the resampling process itself, highlighting the need for targeted refinement in red-band atmospheric processing.

5. Discussion

5.1. Radiometric Calibration Uncertainty Analysis

In this study, radiometric calibration uncertainty is discussed in terms of three primary sources: sensor observation errors, uncertainties associated with the atmospheric correction process, and scale-related effects.

5.1.1. Sensor Radiometric Calibration Uncertainty

The radiometric calibration of Landsat-8 OLI is subject to instrument-level uncertainties, including detector responsivity, electronic noise, optical degradation, and onboard calibration accuracy. The absolute radiometric uncertainty in the visible bands is reported to be approximately 2–3% [42], consistent with values commonly adopted in aquatic remote sensing. This uncertainty is primarily systematic, spatially uniform, and spectrally dependent. In the context of Rayleigh scattering calibration, it sets a baseline for absolute TOA reflectance uncertainty. Since all resampling experiments use the same calibrated Level-1 data, sensor-related errors remain constant across spatial resolutions and methods, mainly affecting the absolute reflectance values while having minimal impact on the relative evaluation of resampling performance.

5.1.2. Atmospheric Correction Uncertainty

In the Rayleigh scattering calibration framework of this study, the uncertainty in the atmospheric correction process mainly arises from two sources: errors in AOD estimation and discretization and interpolation effects in the radiative transfer LUT. Under the marine atmospheric conditions considered here, these two factors are regarded as the dominant contributors to atmospheric correction uncertainty, while other atmospheric parameters are held constant. It should be noted that water vapor and ozone were treated as spatially uniform in this analysis; however, this assumption is reasonable for the selected scenes under marine atmospheric conditions, spatial variability in these constituents may introduce additional uncertainty when the method is applied over broader spatial domains. Assuming the two sources are independent, the total atmospheric correction uncertainty was estimated by combining their individual contributions as:

$$[eqn]$$

AOD Uncertainty

The AOD uncertainty was quantified by introducing a fixed bias (ΔAOD = 0.01) to the AOD input in the 6S radiative transfer model. A maritime aerosol model was selected in accordance with the oceanic calibration environment. With all other atmospheric and geometric parameters held constant, the resulting variation in simulated TOA reflectance yielded the spectral AOD uncertainty for each band.

2.LUT Interpolation Uncertainty

The LUT uncertainty was quantified by comparing TOA reflectance obtained from LUT interpolation with that directly simulated by the 6S model under identical atmospheric and geometric conditions. All LUT input parameters were kept consistent between the two approaches, such that the observed differences primarily reflect LUT discretization and interpolation effects.

As shown in Table 6, the atmospheric correction uncertainty remains below 2% across all spectral bands. LUT interpolation generally contributes slightly more to the total uncertainty than AOD perturbations, particularly at shorter wavelengths (443–486 nm), reflecting the higher sensitivity of shortwave bands to discretization and interpolation effects. In longer wavelengths (561–655 nm), both sources contribute less, resulting in overall smaller uncertainties. These results indicate that the combined atmospheric correction uncertainty is relatively limited across the visible spectrum, supporting the reliability of TOA reflectance retrieval and subsequent water-leaving radiance estimation in the Rayleigh scattering calibration framework.

5.1.3. Spatial Resolution Uncertainty

To evaluate uncertainty related to spatial resolution, the original 30 m images were resampled to coarser resolutions of 500, 1000, and 2000 m using the flux-conserving (FC) method. This approach preserves radiometric flux, so that differences in reflectance mainly result from pixel aggregation and spatial mixing.

For each spectral band, the resampled reflectance was compared with the original 30 m reflectance, which was treated as a reference representation of surface reflectance. The uncertainty was quantified using the relative difference:

$$[eqn]$$

The resulting values were spatially averaged over the study area for each spectral band and resolution. The spatially averaged uncertainties are summarized in Table 7.

Table 7 shows that uncertainty increases as spatial resolution decreases. Shorter wavelengths (443 and 486 nm) exhibit slightly higher uncertainties, whereas longer wavelengths remain more stable. This indicates that image scale is an important factor affecting radiometric measurements at kilometer-level resolutions.

5.2. Flux-Conserving Resampling in Diverse Aquatic Conditions

The superior performance of the flux-conserving (FC) resampling method in optically heterogeneous waters, as established in Section 4, positions it as a promising alternative to interpolation-based techniques. In aquatic remote sensing, however, the efficacy of any spatial processing method is inherently governed by the water body’s inherent optical properties, which vary dramatically across environments. To critically assess the FC method’s general applicability and to address the need for validation under such extreme or distinct optical regimes, we rigorously evaluated its performance across four strategically selected validation sites (GOT-Seaprism, MVCO, Helsinki Lighthouse, and Lucinda). Chosen to represent a broad spectrum of key aquatic optical conditions, these sites provide imagery of distinct, canonical regimes, enabling a targeted assessment of the FC method under different dominant optical processes.

The GOT-Seaprism site in the Gulf of Thailand represents highly turbid nearshore waters with pronounced spatial heterogeneity. The MVCO site along the northeastern coast of the United States corresponds to temperate coastal waters influenced by seasonal phytoplankton growth and variations in suspended particulate matter. The Helsinki Lighthouse site in the Baltic Sea characterizes waters dominated by colored dissolved organic matter (CDOM) absorption. Finally, the Lucinda site in Australian offshore waters represents relatively clear, optically deep waters with a low absorption background, serving as a reasonable proxy for oligotrophic open-ocean optical conditions.

Pixel-level comparisons of water-leaving radiance were performed at three spatial resolutions (500 m, 1000 m, and 2000 m). The RMSE and R^2^ were employed as evaluation metrics to quantify absolute radiometric deviation and the preservation of spatial patterns, respectively. The comprehensive results for each site and resolution are presented in Table 8.

Overall, FC consistently achieves low RMSE and high R^2^ values, even at the coarsest spatial resolution (2000 m). A slight performance degradation occurs as pixel size increases, which reflects the inherent challenge of subpixel heterogeneity. Despite this, the method maintains an R^2^ above 0.83 for all sites. Differences among sites correspond closely to their optical characteristics: high-turbidity nearshore waters (GOT-Seaprism) and CDOM-dominated waters (Helsinki Lighthouse) exhibit relatively higher RMSE, whereas the optically deep and clear offshore site (Lucinda) shows the lowest errors, highlighting the effectiveness of FC resampling under low-absorption, homogeneous conditions.

These results indicate that the FC method reliably preserves radiometric fidelity across diverse water types and spatial resolutions. Its flux-conserving principle mitigates nonlinear radiometric effects induced by pixel aggregation, making it particularly robust in heterogeneous or optically extreme waters. Future work should extend this validation to more dynamic conditions, such as sediment-laden river plumes or intense algal blooms, to further define the boundaries of the method’s robustness.

It is worth noting that the resampling analysis in this study focused on visible bands (443–655 nm), as these wavelengths are most directly relevant to Rayleigh scattering calibration and water-leaving radiance retrieval. Near-Infrared (NIR) and Short-Wave Infrared (SWIR) bands, however, play an equally critical role in atmospheric correction and turbid water detection. Extending this evaluation to include NIR and SWIR bands, therefore, represents an important next step toward assessing the full-spectrum performance of flux-conserving resampling.

Additionally, although scenes with cloud cover below 10% were selected, residual thin clouds may still introduce some uncertainty in radiometric consistency at coarse resolutions, particularly at 2000 m. Future studies could explicitly quantify this effect or adopt advanced atmospheric correction methods, such as neural-network-based approaches near clouds [43], to further improve validation under less ideal conditions.

5.3. Implications for Long-Term Ocean Color Time Series Stability

Building on the multi-scale validation results, the implications for long-term ocean color monitoring are considered here. In addition to its impact on instantaneous radiometric accuracy, the choice of spatial resampling method may also influence the stability of long-term ocean color time series. Many ocean color applications, including water quality assessment, climate-related trend detection, and interannual variability analysis, rely on consistent radiometric behavior across extended temporal records. When multi-year datasets involve varying spatial resolutions or multi-sensor integration, inconsistencies in spatial processing may introduce subtle systematic deviations that propagate into trend analyses.

As demonstrated in Section 4 and Section 5.2, the flux-conserving (FC) method maintains radiometric fidelity across different spatial scales and optical environments. Because it preserves the integrated radiometric flux during spatial aggregation, changes in pixel size primarily reflect spatial mixing effects rather than artificial redistribution of radiance. This property is particularly relevant in long-term monitoring scenarios where spatial resolution may vary between datasets or processing stages.

Interpolation-based methods, while computationally efficient, inherently perform weighted averaging of neighboring pixels. In optically heterogeneous coastal waters, such operations may alter the radiometric balance in a scale-dependent manner. If applied inconsistently across different segments of a temporal record, these effects may introduce artificial discontinuities or scale-dependent variability that is not directly related to environmental change.

In contrast, the FC approach ensures that aggregated reflectance values remain physically consistent representations of the original radiometric signal. Although it does not eliminate subpixel heterogeneity, it minimizes additional methodological uncertainty associated with spatial transformation. This characteristic may contribute to improved comparability across multi-resolution datasets and reduce the risk of bias in long-term trend estimation, thereby enhancing the reliability of climate-quality ocean color data records.

The present study focuses on instantaneous radiometric comparisons. A dedicated multi-temporal evaluation, incorporating extended time series and statistical trend metrics, would be valuable for quantifying the cumulative influence of spatial resampling on long-term ocean color variability. Such analysis would further clarify the role of spatial preprocessing in the interpretation of decadal-scale bio-optical changes and help establish best practices for consistent time series generation.

6. Conclusions

Conventional studies on multi-source aquatic remote sensing have primarily focused on image fusion strategies or bio-optical inversion algorithms, whereas the spatial resampling step preceding fusion has generally been treated as a technical detail. In practice, however, resampling is unavoidable and always occurs before any fusion or retrieval procedure. This study shifts the analytical focus to this often-overlooked preprocessing stage and systematically evaluates how different resampling strategies influence radiometric consistency in water color remote sensing. The main contribution lies in revealing the distortions induced by resampling across spatial scales and in demonstrating the advantages of the flux-conserving method to mitigate these distortions.

The principle of energy conservation is particularly relevant for ocean color remote sensing because aquatic signals are inherently weak and highly sensitive. Water bodies are low-reflectance targets, and radiometric information in the blue and green bands, which is essential for retrieving bio-optical parameters, can be significantly affected even by minor perturbations. Maintaining energy conservation during spatial resampling helps suppress the amplification of such small errors and preserves the fidelity of weak aquatic signals, ensuring more reliable retrieval of water quality indicators.

Quantitative evaluation using Landsat-8 OLI imagery confirms that the FC method consistently outperforms conventional interpolation techniques across all spatial resolutions and spectral bands. Under the most challenging condition examined, namely downsampling from 30 m to 2000 m in the blue band (443 nm), the FC method reduced the mean relative radiometric error to 4.36%, compared with 7.80% for bilinear interpolation, representing an error reduction of approximately 44%. The corresponding standard deviation decreased from 5.31% to 3.59%, indicating substantial improvement in radiometric stability under severe mixed-pixel conditions. Similar advantages were observed in water-leaving radiance retrieval across four geographically diverse AERONET-OC sites, where FC maintained relative errors within ±5% across the blue–green bands for the majority of site and band combinations, and coefficients of determination remained above 0.83, demonstrating robust preservation of spatial patterns and radiometric fidelity even at coarse resolutions.

Multi-scale and multi-regional validation further indicates that the performance of the FC method is largely insensitive to variations in bio-optical conditions. Consistent results were obtained across environments ranging from tropical clear waters to high-latitude turbid systems, highlighting its strong generalization capability and robustness with respect to variations in water constituents and aquatic remote sensing products. The proposed framework also offers a valuable reference for the vicarious calibration of low-resolution sensors and for cross-calibration between sensors with differing spatial resolutions.

Several limitations should be acknowledged when interpreting these results. First, the present analysis is based on Landsat-8 OLI imagery, a sensor characterized by high radiometric stability and a favorable signal-to-noise ratio. Accordingly, the relative performance of the flux-conserving method compared with interpolation-based approaches may vary when applied to sensors with different radiometric characteristics and noise levels. For example, Sentinel-2 MSI, MODIS, and VIIRS differ in signal-to-noise ratio, spatial resolution, and spectral response functions, which may influence the relative magnitude of resampling-induced uncertainty.

Second, although the four AERONET-OC validation sites encompass a broad range of bio-optical conditions—from tropical clear waters to high-latitude turbid systems—they do not fully represent the most optically extreme aquatic environments encountered globally. River-influenced coastal margins with sharp turbidity gradients, glacial meltwater plumes, and nearshore zones with pronounced sub-pixel variability are examples of highly heterogeneous settings in which resampling-induced errors may be amplified.

In addition, the atmospheric correction in this study assumed spatially uniform water vapor and ozone concentrations. While this assumption is reasonable for the selected marine scenes, it may introduce additional uncertainty when the method is applied over broader spatial domains where these constituents exhibit significant spatial variability.

Finally, this study deliberately isolated the effects of spatial resampling under controlled observational geometries. Future work should examine how flux-conserving resampling interacts with viewing geometry effects, including the bidirectional reflectance distribution, and extend validation to more complex atmospheric conditions, such as seasonal cycles, and to additional satellite platforms. Such efforts would provide a more comprehensive assessment of the method’s operational applicability for global aquatic remote sensing.