Sampling Strategies to Improve Passive Optical Remote Sensing of River Bathymetry

Passive optical remote sensing of river bathymetry involves establishing a relation between depth and reflectance that can be applied throughout an image to produce a depth map. Building upon the Optimal Band Ratio Analysis (OBRA) framework, we introduce sampling strategies for constructing calibration data sets that lead to strong relationships between an image-derived quantity and depth across a range of depths. Progressively excluding observations that exceed a series of cutoff depths from the calibration process improved the accuracy of depth estimates and allowed the maximum detectable depth (dmax) to be inferred directly from an image. Depth retrieval in two distinct rivers also was enhanced by a stratified version of OBRA that partitions field measurements into a series of depth bins to avoid biases associated with under-representation of shallow areas in typical field data sets. In the shallower, clearer of the two rivers, including the deepest field observations in the calibration data set did not compromise depth retrieval accuracy, suggesting that dmax was not exceeded and the reach could be mapped without gaps. Conversely, in the deeper and more turbid stream, progressive truncation of input depths yielded a plausible estimate of dmax consistent with theoretical calculations based on field measurements of light attenuation by the water column. This result implied that the entire channel, including pools, could not be mapped remotely. However, truncation improved the accuracy of depth estimates in areas shallower than dmax, which comprise the majority of the channel and are of primary interest for many habitat-oriented applications.


Introduction
Beginning in the 1990s, remote sensing emerged as a powerful means of characterizing river channel morphology and in-stream habitat [1][2][3][4][5][6]. Buoyed by advances in sensors and software, Marcus and Fonstad [7] and Carbonneau et al. [8] highlighted the potential for remote sensing to facilitate river research and management. Whereas earlier methodological studies focused on algorithm development and testing at the reach scale [9][10][11][12], more recent work is broader in scope, encompassing long river segments or entire watersheds [13][14][15]. This kind of large-scale remote sensing could support a broad range of applications, such as assessing geomorphic impacts of extreme flow events [16] and quantifying habitat availability for critical species [17]. In a hydrologic context, measuring river discharge from remotely sensed data is an important research objective and several approaches have been explored: hydraulic approximations based on observations of width, slope, and flow resistance [18], estimating surface flow velocities and channel bathymetry from optical [19] where R(λ 1 ) and R(λ 1 ) are reflectances, radiances, or raw digital numbers recorded in numerator and denominator bands centered at wavelengths λ 1 and λ 2 , respectively. The Optimal Band Ratio Analysis (OBRA) algorithm involves performing regressions of X vs. d for all possible band combinations and identifying the (λ 1 , λ 2 ) pair that yields the highest coefficient of determination R 2 . Legleiter et al. [29] showed that OBRA provides reliable depth estimates in the presence of variable bottom types, optical properties, and surface reflectances, but also cautioned that spectrally based bathymetric mapping is only feasible under certain conditions, with the primary requirement being clear, relatively shallow water. Subsequent work indicated that incorporating an X 2 term could improve depth estimates in deeper channels [4,30]. Although OBRA is designed to exploit the detailed spectral information provided by hyperspectral images or field spectra, Legleiter et al. [29] and Legleiter and Overstreet [4] showed that accurate maps of river bathymetry also can be produced by applying the OBRA algorithm to multispectral satellite images with fewer, much broader bands. Similarly, depths can be estimated using simple log-transformed band ratios calculated from standard, publicly available red, green, blue, and possibly near-infrared (RGB+NIR) digital aerial photography such as that available through the National Agricultural Imagery Program (NAIP; [31]).
OBRA has become an established method of bathymetric mapping, but the technique is empirical and requires ground-based depth measurements for calibration. This reliance on field data is an important limitation that also raises the question of how calibration data influence depth retrieval. Biases introduced during field data collection could propagate through OBRA and affect image-derived depth estimates. For example, real-time kinematic (RTK) GPS-equipped wading surveys along shallow channel margins are often coupled with acoustic Doppler current profiler (ADCP) or echo sounder measurements from deeper areas. Because manual surveys provide fewer measurements per unit time than boat-based instruments, shallow depths might not be included in the merged calibration data set in proportion to their actual areal extent and thus could be under-represented. An X vs. d relation calibrated by using field data that do not represent all depths present in the river in the proper proportions might not provide reliable estimates across the full range of depths.
The effects of biased sampling on depth retrieval have not been examined quantitatively, but our experience on a number of bathymetric mapping projects [32], as well as statistical reasoning, suggests that the distribution of depths used for calibration influences the X vs. d relationship.
In addition, the statistical approach used to fit such a relation and the functional form used to describe the relationship are important aspects of the overall calibration process. In this study, we focus on sampling while holding the curve-fitting method (regression) and type of relationship (quadratic) constant. Although both the independent (X) and dependent (d) variables are subject to uncertainty, regression is an appropriate method of model fitting because the coefficients are not of direct interest and the goal is to predict d from X [33,34]. The X vs. d relation could be represented via other, nonlinear functional forms, such as exponential or power law models, but herein we use the same kind of second-order polynomials developed for estimating depth in previous studies [4,30]. We concentrate on sampling as an initial step toward a better understanding of the calibration process because the way in which a calibration data set is assembled will influence all subsequent stages of that process. Examining more advanced methods of curve-fitting and considering alternative functional forms, given a calibration data set, are important topics for future research.
One consequence of sampling, particularly for OBRA relations that include an X 2 term such that d is not monotonically related to X, is that a lack of shallow field observations can lead to minimum predicted depths greater than the smallest actual depths; the resulting bathymetric maps are too deep along the edges of the channel. This scenario is illustrated in Figure 1, an example from the Deschutes River in Oregon, USA. The quadratic X vs. d relation shown in Figure 1a reaches a minimum at a depth of 0.73 m, with lesser values of X resulting in deeper rather than shallower depth estimates. Applying this OBRA equation to the hyperspectral image thus leads to a bathymetric map with unrealistically deep water near the banks (Figure 1b). The observed vs. predicted regression used for validation highlights the absence of image-derived depth estimates shallower than 0.73 m despite numerous field observations in this range (Figure 1c).
In addition to the issues associated with shallow channel margins, estimating depth in pools also can be problematic. Whereas inaccurate depth estimates in shallow areas are most likely sampling-related artifacts of calibration, systematic underestimation of pool depths is a consequence of the physics governing the interaction of light and water: beyond a certain maximum detectable depth, denoted by d max (λ), further increases in depth no longer produce a change in water-leaving radiance of sufficient magnitude to be resolved by the imaging system [35]. Both environment and instrumentation thus constrain spectrally based bathymetric mapping, with d max (λ) depending not only on the optical characteristics of the river but also the technical specifications of the sensor.
Philpot [35] expressed these concepts in terms of the ratio of the noise-equivalent delta radiance ∆L N (λ), an index of detector sensitivity, to the radiance reflected from the bottom, L B (λ). In clear-flowing streams and for wavelengths dominated by pure water absorption, the amount of radiance propagating to the streambed and reflecting back upward decreases as depth increases. For a given sensor, (i.e., value of ∆L N (λ)) the ratio ∆L N (λ)/L B (λ) thus increases with increasing depth. The detection limit occurs when the difference between the total at-sensor radiance and the radiance from a hypothetical infinitely deep water column is equivalent to the sensor's ∆L N (λ), implying that the streambed is essentially no longer visible.
Legleiter and Roberts [36] introduced a forward image modeling framework for simulating these effects, but ideally such a model would be parameterized based on a priori knowledge of river bathymetry, bottom reflectance, and water column optical properties. Specialized equipment is required to obtain this kind of optical field data, however, and performing a detailed channel survey to enable forward image modeling might obviate the need for remote sensing. Similarly, although measurements of downwelling spectral irradiance E D (λ) at different depths within the water column can be used to calculate d max (λ) [11,37], such data are rarely available. An alternative approach to identifying d max (λ) directly from an image thus would be valuable and could help to avoid spurious estimates of pool depths and derived quantities such as erosion volumes. The maximum detectable depth also is pertinent to sampling and calibration because any field-based depth measurements in excess of d max (λ) will not help to inform the X vs. d relation but rather diminish the strength of the relationship. Knowledge of d max (λ) thus could be used to constrain the range of depths that should be included in the calibration process. Figure 1. Example of Optimal Band Ratio Analysis (OBRA)-based depth retrieval from the Deschutes River using a quadratic relation between the image-derived quantity X and flow depth d. (a) calibration plot showing the non-monotonic relationship between X and d due to the X 2 term. The X vs. d equation represented by the red line reaches a minimum at d = 0.73 m; (b) depth map produced from a hyperspectral image based on the quadratic X vs. d relation shown in (a). Note the depths along the channel margins do not go to to zero but rather the minimum estimated depth of 0.73 m. Flow is from left to right; (c) observed vs. predicted (OP) regression used to validate depth estimates, showing the absence of depth estimates less than 0.73 m. The best-fit regression line is shown in red. d f refers to the field-surveyed depth, d i to the image-derived depth estimate, n is the sample size, and σ R is the standard error of the regression in units of m.
Given the demonstrated potential for, and increased interest in, remote sensing of rivers, our overall goals are to develop and test refined methods for estimating water depth from passive optical image data and to better understand the inherent limitations of these techniques. More specifically, this study is motivated by a question that has not been considered previously: how are image-derived depth estimates affected by the distribution of field measurements used for calibration? Using field spectra and hyperspectral image data from two distinct river systems, we pursue the following research objectives:

1.
Introduce sampling strategies for improving the accuracy of depth estimates by (a) truncating the calibration data set to exclude field observations beyond a specified maximum cutoff depth; and/or (b) stratifying the calibration data set such that all ranges of depths present in the river of interest are represented in equal proportions.

2.
Demonstrate how this type of sampling can be used to estimate the maximum detectable depth directly from an image, without field measurements of water column optical properties.

3.
Compare empirical, image-based estimates of d max to theoretical calculations based on irradiance profiles recorded in the field.

Study Areas and Field Data Collection
We developed and tested methods for spectrally based depth retrieval using field measurements and hyperspectral images from two distinct river systems: the Deschutes River in Oregon, USA, and the Snake River in Wyoming, USA (Figure 2a). These two rivers differ from one another in terms of their morphology, clarity, and available data and thus pose disparate challenges for bathymetric mapping. The Deschutes is narrower, deeper, and less clear than the Snake, implying that estimating depths could be more problematic on the Deschutes. However, we obtained an essentially continuous, field-based bathymetric survey of the Deschutes to calibrate image-derived depth estimates. For the Snake, in contrast, only a more typical, relatively sparse field data set was available for calibration. All of the data used in this study are included in data releases available through the USGS ScienceBase catalog [38,39].
On the Deschutes, we focused on a 3.67 km reach near Bend, between Benham Falls and Dillon Falls, where the channel is approximately 45 m wide, with a gentle gradient, low velocities, and depths up to 8.87 m ( Table 1). The river's course is dictated by lava flows and basalt is exposed on the bed and banks in several locations. The impetus for work on the Deschutes was the need to characterize in-stream habitat for the Oregon Spotted Frog (Rana pretiosa) to help guide flow management conducive to the species' well-being [40]. We surveyed the channel using an ITER Systems Bathyswath-1 468 kHz multi-beam echo sounder (Saint-Jorioz, France) deployed from a cataraft. By making several parallel up-and downstream passes, we obtained essentially continuous, bank-to-bank bathymetric coverage ( Figure 2b). Post-processing of the raw survey data involved de-spike and smoothing filters and interpolation onto a 0.5 m grid; for further details, see the metadata associated with [38]. This exhaustive field data set provided information on the actual distribution of depths within the channel and thus was ideal for examining the effects of sampling strategy on spectrally based depth retrieval. Although the Deschutes did not convey high concentrations of suspended sediment during our study, organic materials reduced water clarity and the channel appeared visually dark, with a low-reflectance substrate composed of basaltic sands and silts.
On the Snake River, we examined a 1.2 km reach in Grand Teton National Park with a mean width of 56 m. This site has been the subject of several channel change studies [41][42][43] and a testing ground for various methods of mapping depth [4,24,25,31], bed elevation [44], and bottom composition [45] from remotely sensed data. We obtained depth measurements through a combination of wading surveys along channel margins and bar surfaces and ADCP cross-sections collected from a kayak for deeper areas (Figure 2c). The wading measurements were concentrated in shallow portions of the channel and are evident in Figure 2c as isolated, more widely spaced, distinct points near the banks that are offset from the ADCP transects that spanned most of the channel width. While wading, we measured depths directly on a top-setting rod. ADCP depth measurements were obtained using a SonTek RiverSurveyor M9 (San Diego, CA, USA) configured to record the vertical beam as the primary depth reference; for further details, refer to the metadata accompanying [39]. The Snake had a mean depth of 0.98 m (Table 1), shallower than the Deschutes, and was clearer, making the former river more favorable for depth retrieval. The hybrid wading/ADCP data set, however, implied that spectrally based bathymetric mapping could be compromised by sampling-related issues, particularly under-representation of shallow depths (Table 1).  In addition to depth measurements, we collected field spectra from both rivers using an Analytical Spectral Devices FieldSpec3 (Westborough, MA, USA) spectroradiometer mounted on a raft and operated in reflectance mode; see [45] for further details on this instrument and our protocol. The depth at each spectral measurement location was interpolated from data collected simultaneously by a single-beam echo sounder. These field spectra were essentially continuous, with a 1 nm sampling interval from 400 to 900 nm, and thus provided greater spectral resolution than the hyperspectral images. Moreover, using field spectra allowed us to examine the effects of sampling strategy on depth retrieval in the absence of geo-referencing issues, atmospheric effects, and other factors that complicate the use of image data.
We also used the spectroradiometer to measure the amount of downwelling spectral irradiance E D (λ) propagating to different depths within the water column. To obtain these irradiance profiles, we mounted an upward-facing fore-optic with a cosine receptor on the end of a measurement crane that allowed us to efficiently raise and lower the instrument. A pressure transducer was placed on the crane next to the fore-optic to determine the depth of each measurement. We used the irradiance spectra to calculate the diffuse attenuation coefficient K d (λ) following Mishra et al. [37], which in turn was used to estimate the maximum detectable depth for a specified ∆L N (λ)/L B (λ) ratio following Philpot [35]. This analysis resulted in field-based estimates of d max (λ) for comparison with d max values inferred directly from image data.

Remotely Sensed Data and Image Processing
Hyperspectral image data were acquired from both the Deschutes and Snake Rivers by the ITRES Compact Airborne Spectrographic Imager (CASI) 1500H (Calgary, Alberta, Canada), a pushbroom scanner with 1500 across-track imaging pixels and up to 288 spectral bands [46]. Flight planning for this type of sensor involves a compromise between spatial and spectral resolution [44], and, in this study, we obtained images with 0.5 m pixels and 48 bands equally spaced between 368 and 1039 nm. For the Deschutes, we discarded the shortest and longest wavelength bands due to excessive noise at both ends of the spectrum, which is common in hyperspectral images of rivers due to atmospheric effects for λ < 400 nm and strong absorption by pure water for λ > 1000 nm, and retained 42 bands ( Table 1). The raw CASI data were radiometrically corrected to obtain spectral radiance units of µW cm −2 sr −2 nm −1 using software and calibration coefficients provided by the manufacturer. We used the resulting radiance image for the Deschutes, but for the Snake we used a reflectance image produced via the ATCOR atmospheric correction program (version 4, ReSe Applications LLC, Wil, Switzerland) [47]; for further details, refer to [44] and the metadata for [39]. Both images were geo-referenced based on trajectories recorded by a GPS and inertial motion unit on the aircraft during each flight. The resulting images were accurately co-registered with surveyed ground control targets and no spatial adjustments were required.
For both rivers, the image processing workflow involved defining a channel mask by applying a NIR band threshold and manually editing the initial mask as needed, clipping the field data to this mask, and applying a 3 by 3-pixel Wiener spatial smoothing filter to the in-stream portion of each image. For the Snake River, pixel-scale mean depths were calculated by assigning each field measurement point to the closest pixel center and then averaging all of the depths assigned to a given pixel. For the Deschutes, this averaging step was unnecessary because the multi-beam echo sounder data were gridded to the same 0.5 m resolution as the CASI image data. We estimated depths by performing OBRA based on calibration data sets obtained via the sampling strategies described below.

Sampling Strategies for Improved Depth Retrieval
The spectrally based bathymetric mapping techniques outlined herein are based upon the following premise: because the nature and strength of a relationship between an image-derived quantity and depth largely depends on the data used to calibrate that relation, field data should be sampled to optimize the calibration process. Even if the field effort yielded a data set that included measurements greater than the maximum detectable depth and/or favored certain ranges of depths, we reasoned that potential biases might be avoided by selectively sampling field observations to obtain a calibration data set conducive to accurate depth retrieval across the full range of depths. Essentially, our approach involves extracting a subset of the available field data and then using this subset to calibrate a relationship between depth and reflectance via OBRA. This paper introduces two methods for implementing this workflow and evaluates these sampling strategies using field spectra and hyperspectral images from the Deschutes and Snake Rivers.

Optimal Band Ratio Analysis (OBRA) of Progressively Truncated Input Depths (OPTID)
This sampling strategy is intended to enhance depth retrieval by excluding field observations that exceed the maximum depth detectable by a particular sensor. OPTID thus provides a means of inferring d max directly from an image, without field measurements of water column optical properties. OPTID involves removing field observations greater than a specified cutoff depth from the data set used to calibrate an X vs. d relation via OBRA. The procedure is iterative, repeated for each one of a series of different cutoff depths, with a separate OBRA conducted at each step. In this study, we progressively truncated the input depths in 5 cm increments, starting with the maximum depth in each of the original field data sets and continuing down to approximately 0.5 m. We quantified the effects of this sampling strategy in terms of both the R 2 values of X vs. d relations calibrated via OBRA and the R 2 values of observed vs. predicted (OP) regressions used to validate image-derived depth estimates [48]. The latter were based on a 50% random sample of the original field data withheld from the calibration process.

Stratified Optimal Band Ratio Analysis (SOBRA)
SOBRA is a strategy for sampling a collection of field measurements to obtain a calibration subset in which all ranges of depths are represented equally, even if the original distribution is far from uniform. This procedure thus provides a means of accounting for field methods that tend to preferentially, or at least more efficiently, sample certain areas of the channel, possibly biasing the calibration process and resulting image-derived depth estimates. By separating the actual spatial frequency distribution of depths from the purely statistical distribution of data used for calibration, SOBRA prioritizes establishing an X vs. d relationship that is applicable across the full range of depths, regardless of their abundance within the channel. This approach is intended to produce an X vs. d relation well-suited for application to shallow, intermediate, and deep areas rather than an equation that performs well for portions of the river but yields less reliable depth estimates along channel margins and/or in pools.
SOBRA involves binning the depth measurements and then retaining the same number of samples from each bin. Input parameters for SOBRA include the number of depth bins and the lower limit of the deepest bin, specified as a percentile of the distribution of depths in the original field data set. Histograms for the original, binned, and final SOBRA calibration data sets for the Deschutes and Snake rivers are illustrated in Figure 3. We used 10 bins, with the lower bin limits equally spaced between the minimum depth in the field data set and the 95th percentile, denoted by d 95 ; all depths greater than d 95 are assigned to the last (deepest) bin (Figure 3b,e). The bin containing the smallest number of measurements is used to establish the sample size and the same number of observations are randomly selected from each of the other bins (Figure 3c,f). Once established, the stratified calibration data set is used as input to OBRA. OBRA and OP regression R 2 values were used to quantify the effects of stratification on depth retrieval accuracy.

Estimating the Maximum Detectable Depth
We used OPTID results to infer the maximum detectable depth directly from hyperspectral image data and field-based reflectance spectra from the Deschutes and Snake rivers. For each cutoff depth, we performed observed vs. predicted (OP) regressions based on field data reserved for validation and plotted the resulting R 2 values, along with those from the OBRA regressions used to calibrate X vs. d relationships, against the greatest depth retained in the calibration data set. We hypothesized that beyond a certain cutoff depth, the OP and OBRA R 2 values would begin to decrease as field observations greater than d max were included in the data set. An inflection point on a plot of R 2 vs. cutoff depth thus could provide an indication of d max . Moreover, the value of X corresponding to d max could be calculated by rearranging the X vs. d relation from the OBRA for that cutoff depth and then used as a threshold to identify image pixels deeper than d max .
To evaluate this new image-driven approach, we compared OPTID-based estimates of d max to values derived from irradiance profiles measured in the field, calculated following Philpot [35] as For consistency with earlier work [4], we used ∆L N (λ)/L B (λ) values of 0.1 and 0.01 to illustrate the effects of a greater bottom contrast between the streambed and water column and/or a more sensitive detector, either or both of which would correspond to smaller values of ∆L N (λ)/L B (λ) and thus greater d max (λ). Whereas d max (λ) was calculated from the irradiance profiles for all wavelengths, OPTID-based estimates of d max were based on the band combination identified by OBRA and did not vary with wavelength. To compare the two approaches, we thus used d max (λ) values calculated from the irradiance profiles for the numerator and denominator wavelengths identified by OBRA for the cutoff depth with the highest OBRA R 2 .

OPTID of Hyperspectral Image Data and Field Spectra
OBRA of Progressively Truncated Input Depths (OPTID) seeks to improve depth retrieval by removing field measurements that exceed each one of a series of cutoff depths from the calibration data set. The first panel of Figure 4 shows the OBRA matrix for Deschutes image spectra based on the original field data, followed by OBRA matrices for OPTID cutoff depths from 0.57 m, for which only 9.8% of the field data were included, up to 8.57 m, for which 99.9% of the original depth measurements were retained. The broader areas of warm yellow, orange, and red tones across the top row of Figure 4 indicate that the strength of the optimal X vs. d relation increased as the cutoff depth increased from 0.57 to 3.57 m. At greater cutoff depths, however, including deeper field observations in the calibration data set reduced the OBRA R 2 . Because nearly all of the field data were retained for the largest cutoff depth of 8.57 m, the OBRA R 2 decreased to a value similar to that for the original data set. The OBRA matrices for all 168 cutoff depths applied to the Deschutes River are compiled in an animation provided as Supplementary Information online. Video S1 more dynamically illustrates how the strength of the X vs. d relation (represented by higher OBRA R 2 values and warmer tones in the matrix plots) initially increased with cutoff depth, reached a plateau, and then decreased for cutoff depths greater than 3.57 m. Figure 5 depicts a similar pattern for shallow cutoff depths on the Snake River, with OBRA R 2 increasing and the range of useful wavelength combinations expanding as cutoff depth increased up to about 2 m. OBRA R 2 then stayed relatively constant for cutoff depths approaching 3 m, rather than decreasing for larger cutoff depths as occurred on the Deschutes. Even when all of the field data were included (i.e., the original data set), the first panel in Figure 5 indicates that the relationship between X and d remained strong (R 2 = 0.79) across a broad range of visible band ratios. OBRA matrices for all 51 cutoff depths used on the Snake River are combined in Video S2, which shows the strength of the X vs. d initially increasing rapidly with cutoff depth, stabilizing around 2 m, and then remaining steady as the cutoff depth approached the largest observed depth. Sequences of OPTID matrices like those shown in Figures 4 and 5 also illustrate how the range of useful wavelengths changes as the range of depths included in the calibration data set changes. The optimal numerator and denominator bands varied with depth due to spectral differences in the absorption of light by pure water, which is relatively weak at blue and green wavelengths but increases rapidly in the red and NIR. Consequently, as the cutoff depth applied to the Deschutes data set decreased from 8.57 to 3.57 m, OBRA R 2 values for NIR denominator bands up to 800 nm increased for green numerator bands from 500 to 600 nm. The wavelength combinations identified by OBRA for the various cutoff depths used at each site are shown in Figure 6. For the Deschutes, λ 1 and λ 2 remained steady at a green band (520 nm) and a band on the edge of the NIR (700 nm) for relatively shallow cutoff depths up to 4 m. For cutoff depths greater than 4 m, strong absorption of NIR light by the deeper water caused shorter, less strongly absorbing red denominator wavelengths closer to 620 nm to be selected along with a longer numerator wavelength. On the Snake River, the optimal band ratio consisted of two adjacent bands near 600 nm and did not vary for cutoff depths greater than 0.6 m. To summarize how depth retrieval accuracy varied with OPTID cutoff depth, we plotted OBRA R 2 and observed vs. predicted (OP) regression R 2 against cutoff depth for both field sites and data types. Figure 7a,b show that for OPTID of image spectra from the Deschutes and Snake rivers, OBRA R 2 initially increased with cutoff depth as a greater number of field measurements spanning a broader range of depths were retained for calibration. For the Deschutes, OBRA R 2 reached a maximum at a cutoff depth of 3.57 m and then decreased for greater depths. The Snake exhibited a different trend: OBRA R 2 reached a plateau at approximately 1.8 m but did not decrease for greater cutoff depths. Unlike the image spectra, OBRA R 2 for OPTID of field spectra initially decreased with cutoff depth before increasing when more, deeper observations were included in the calibration data sets (Figure 7c,d). OP regression R 2 was generally lower than OBRA R 2 and also varied less as a function of cutoff depth than OBRA R 2 because the bands selected by OBRA were stable across a range of cutoff depths and thus the predicted depths at validation points varied little if at all.

Original calibration data
We assessed the sensitivity of the OPTID algorithm to sampling strategy and sample size as described in the online Supplementary Information (Text S1). Plots of OBRA R 2 vs. cutoff depth followed very similar trends for both stratified and non-stratified samples and for sample sizes that varied as a function of cutoff depth or were held constant for all cutoff depths ( Figure S1). These results imply that OPTID is robust to sampling effects and can be applied directly to the original data set, without any kind of special sampling required.

SOBRA of Hyperspectral Image Data and Field Spectra
To provide a baseline for comparison with SOBRA, we performed standard (non-stratified) OBRA for CASI images and field spectra from both study sites. An example of the output from standard OBRA of the Deschutes image is shown in Figure 1, which indicates moderate OBRA calibration R 2 and OP validation R 2 values of 0.59. Moreover, the quadratic OBRA relation in Figure 1a led to unrealistically high depth estimates along channel margins. Although using a different functional form to describe the X vs. d relation could mitigate this issue, a possibility that should be evaluated in future studies, a stronger relationship between reflectance and depth also can be obtained via stratified sampling. For example, generating a stratified calibration data set prior to OBRA of the Deschutes River image yielded a modest improvement in predictive power, with an OBRA R 2 of 0.64 for 10 depth bins ( Table 2). The results of similar analyses for the Deschutes field spectra and both image and field spectra from the Snake are provided online as Supplementary Information. For the Deschutes field spectra, stratification produced an increase in OBRA R 2 , from 0.79 to 0.86 (Table S1). Differences between standard and stratified OBRA (SOBRA) were less pronounced for the Snake, but stratification improved OBRA R 2 from 0.79 to 0.80 for image data and from 0.85 to 0.88 for field spectra (Tables S2 and S3).

Image-Based Inference of the Maximum Detectable Depth
OPTID could provide a means of estimating the maximum detectable depth directly from an image as an inflection point on a plot of OBRA R 2 vs. cutoff depth. To assess the feasibility of this approach, we compared d max values inferred via OPTID of image spectra to d max (λ) calculated from irradiance profiles measured in the field ( Of the four data sets we considered, only the Deschutes image spectra yielded a plot of OBRA R 2 vs. cutoff depth with a clear inflection point that might provide an indication of the greatest depth the CASI sensor could detect (Figure 7a). For the Snake image and field spectra from both sites, no such maximum was evident and d max could not be inferred unambiguously (Figure 7b-d). For the Deschutes image, d max (λ 1 ) for the 519 nm band selected as the optimal numerator was 2.78 m for an assumed ∆L N (λ)/L B (λ) of 0.1, comparable to but shallower than that estimated from OPTID of the image. For a more sensitive instrument with ∆L N (λ)/L B (λ) = 0.01, d max (λ 1 ) would be 5.55 m. For field spectra from the Deschutes, d max (λ 1 ) was 3.56 m for ∆L N (λ)/L B (λ) = 0.1, very close to the image OPTID-based estimate, and 7.12 m for ∆L N (λ)/L B (λ) = 0.01. For the clearer-flowing Snake River, d max could not be estimated from OPTID of the image, but calculations based on measured irradiance profiles suggested maximum detectable depths of 8.99 and 17.98 m for ∆L N (λ)/L B (λ) values of 0.1 and 0.01, respectively, at λ 1 = 583 nm, the numerator band selected by OBRA for both image and field spectra.

OPTID: A Trade-Off between Signal Strength and Saturation
The sampling strategies introduced herein are designed to facilitate spectrally based mapping of river bathymetry from passive optical image data. Optimal Band Ratio Analysis (OBRA) of Progressively Truncated Input Depths (OPTID) also provided insight regarding the range of depths over which image-derived depth estimates might be reliable. For example, the OBRA matrices based on the Deschutes image shown in Figure 4 indicated that excluding field observations deeper than 3.57 m from the calibration data set led to X vs. d relations with greater predictive power. The decrease in bathymetric accuracy that occurred when measurements with d > 3.57 m were retained suggests that the CASI sensor could not detect areas deeper than about 3.5 m on this river. The inclusion of deeper measurements adversely affected the calibration process because increases in depth beyond 3.5 m were not associated with changes in radiance resolvable by the sensor and thus weakened the relationship between X and d.
Conversely, for both the Deschutes and Snake Rivers, if the range of depth measurements retained for calibration was too small, OBRA R 2 values were very low (Figures 4 and 5). This reduction in predictive power for calibration data sets consisting exclusively of shallow observations suggests that the X vs. d relationship deteriorated as the range of depths used for calibration decreased. We attribute this trend to a decline in the signal to noise ratio for OBRA, where the signal of interest is variation in radiance as a function of depth. Noise arises from the inherent variability associated with differences not only in depth but also bottom composition, surface reflectance, and water column optical properties. In addition, various environmental and instrumental effects can introduce noise to image spectra. If the range of d is too small, the corresponding radiance signal will not stand out above this background noise and any relationship between X and d will be weak; an equation based exclusively on shallow depths might not perform well across the full range of depths. This concept is depicted graphically in the second panel of both Figures 4 and 5, where the OBRA matrices consist entirely of dark, cool tones representing weak X vs. d relations based on only a small range of shallow depths.
In general, a compromise must be reached between excluding depths beyond d max and preserving enough variation in d to provide an X vs. d signal strong enough to overcome the noise that would overwhelm any relationship between depth and reflectance if the calibration data set does not include a sufficient range of depths. The OPTID summary in Figure 7 illustrates this trade-off between signal strength and saturation. At the lowest cutoff depths, all but the shallowest field observations were removed from the calibration data sets. For such a small range of depth, variations in reflectance driven by changes in depth were subtle and obscured by changes in reflectance associated with bottom type, sun glint, optical properties, and other factors. As cutoff depth increased, OBRA R 2 increased as more variation in d provided a stronger signal that drowned out this background noise and allowed more robust X vs. d relations to be established. For the Deschutes, OBRA R 2 was maximized for a cutoff depth of 3.57 m but decreased for greater depths as the calibration data set began to include a larger number of depths beyond what the CASI imaging system could detect. For the Snake, in contrast, OBRA R 2 stopped increasing beyond a cutoff depth of 1.8 m but remained steady and did not decrease for greater cutoff depths. This result implies that the CASI sensor's d max was not exceeded on the Snake due to the greater water clarity and shallower depths in this river than on the Deschutes. In addition to turbidity, two other factors contributed to the ability of the CASI sensor to detect the full range of depths present on the Snake but not on the Deschutes. First, the Deschutes channel included areas that were much deeper than any pools within our study reach on the Snake River and might not have been fully mappable even if the water in the Deschutes were clearer. Second, the Deschutes is within a region of volcanic lithology and flows directly over dark basalt in some places, whereas the Snake has a more typical, brighter gravel-cobble bed. The lower bottom reflectance of the Deschutes thus might have reduced the d max for this river as well.
The OPTID of field spectra summarized in Figure 7c,d provided additional evidence of the importance of having a sufficiently strong signal to establish a relation between depth and reflectance. In contrast to the image spectra, for field spectra OBRA R 2 initially decreased with cutoff depth up to 1.37 m on the Deschutes and 0.70 m on the Snake before increasing when more, deeper observations were included in the calibration data sets. When the range of depths was small, the reflectance signal associated with variations in d was weak and obscured by the inherent variability of field spectra from natural channels. Only when a broader range of depths was considered did stronger X vs. d relations emerge. In addition, the number of observations available was much less for field spectra than for image spectra, especially on the Deschutes, and the low R 2 values also were a consequence of small sample sizes when only shallow field spectra were retained.
For image spectra from both rivers, the OP regression R 2 used as a validation metric was less than OBRA R 2 for all but the shallowest cutoff depths due to the presence of more, greater depths in the validation data set, many of which might have exceeded d max . Unlike OBRA R 2 , OP R 2 did not decrease when greater depths were included in the calibration data set. Whereas the OBRA performed for each OPTID increment was based on a calibration sample that was screened to exclude depths exceeding a certain cutoff, OP regression was based on 50% of the original field data, sampled at random and not screened. As a result, OP R 2 increased as deeper field observations were retained and the calibration data set became more similar to the unscreened validation data set. An increased number of deep field measurements in the calibration data set caused the X vs. d relation to become more appropriate for the greater depths that were more common in the validation data set.
Another factor contributing to the stabilization of OP R 2 values is summarized in Figure 6, which indicates that the numerator and denominator wavelengths identified as optimal (i.e., the band ratios yielding the strongest relationships between depth and reflectance) remained stable across a range of OPTID cutoff depths, provided a sufficient range of depths was retained. Because the same pair of bands was used, the predicted depths at the validation points changed very little and OP R 2 remained steady. For portions of the Deschutes less than 4 m deep, depth information primarily came from the denominator band at 700 nm due to the sensitivity of near-infrared (NIR) reflectance to variations in depth in shallow water. For OPTID cutoff depths greater than 4 m, however, the NIR radiance signal began to saturate and a shorter, red wavelength was selected for λ 2 . This denominator band was not as responsive to changes in depth in shallow water but allowed for greater penetration. For the Snake, which was less deep and clearer than the Deschutes, OBRA identified the same two bands, both near 600 nm, for all cutoff depths greater than 0.6 m. These results suggest that the band ratios and X vs. d regression coefficients identified via OBRA were stable across a range of cutoff depths.

Stratified Optimal Band Ratio Analysis (SOBRA)
Relative to standard OBRA of the original field data, OBRA R 2 increased slightly when the calibration data were stratified into a set of depth bins that all contained the same number of observations. This improvement occurred for both the Deschutes and Snake rivers and for both image data and field spectra, suggesting that in some cases depth retrieval can be enhanced by stratifying to obtain calibration data sets that represent all depth ranges in equal proportions.
Whereas OPTID is an iterative procedure that involves conducting a separate OBRA for each one of a sequence of cutoff depths, SOBRA first establishes a stratified calibration data set and then performs a single OBRA. The SOBRA algorithm is thus more computationally efficient than OPTID, which could be an important consideration for data sets consisting of many spectral bands and/or a large number of field-based depth measurements. SOBRA also is well-suited for application to hybrid field data sets, such as the wading and ADCP surveys from the Snake River, because stratification accounts for differences in point density and ensures that all depth ranges are represented equally. OPTID is more appropriate for field data acquired by a single method, such as the multibeam echo sounder survey of the Deschutes, and/or when inferring the greatest remotely sensible depth is a primary objective. Conversely, SOBRA might yield more realistic depth estimates along shallow channel margins by ensuring that small depths are not under-represented in the calibration data set. Stratified sampling could help to avoid the over-prediction of shallow depths that can occur when a quadratic X vs. d relation is disproportionately influenced by field data from deeper areas of the channel. OPTID and SOBRA also could be used together, with the d max inferred via OPTID used to truncate the data set prior to stratification for SOBRA. This approach would provide a highly robust calibration data set by not only excluding depths greater than the sensor's d max but also ensuring that shallow field observations are well-represented. An important implication of this study is that successful depth retrieval across the full range of depths requires at least some calibration data from shallow areas of the channel, which might involve multiple methods of field data collection, such as the hybrid wading/ADCP survey we performed on the Snake River.

Alternative Approaches to Estimating the Maximum Detectable Depth
Knowledge of the maximum depth a particular imaging system can detect in a specific stream can help to avoid spurious estimates of pool volumes, scour depths, and other geomorphic and ecological parameters. For this reason, a direct, image-based method of estimating d max would be valuable. The OPTID algorithm described herein can provide insight as to which depths might be outside a sensor's dynamic range. A plot of OBRA R 2 vs. the maximum depth retained in the calibration data set (Figure 7) can be used to identify the depth beyond which the strength of the relationship between X and d begins to deteriorate due to saturation of the radiance signal in deeper water.
In this study, only one of the four data sets we examined, image spectra from the Deschutes River, exhibited an obvious inflection point from which d max could be inferred. For field spectra from the Deschutes, OBRA R 2 reached a plateau at about 3.7 m, comparable to the 3.57 m at which OBRA R 2 began to decrease for the image spectra, but did not decrease significantly for depths up to 5.37 m. This result could imply that the field-based spectroradiometer was more sensitive than the airborne imaging system and/or that data collected from a raft were less subject to environmental noise (e.g., atmospheric effects) that would increase a sensor's ∆L N (λ) and thus reduce its sensitivity. Another factor that might have contributed to the discrepancy between OPTID results for image vs. field spectra could be the greater abundance of deep field observations in the multibeam echo sounder survey used to analyze the image. Analysis of the field spectra, in contrast, was based on depths from a single beam echo sounder that provided only point measurements, not exhaustive coverage of the channel like the multibeam system. In any case, the analysis summarized in Figure 7a indicates that for this image of the Deschutes, depth estimates greater than 3.57 m should be treated with caution or disregarded altogether. Because our multibeam survey indicated that approximately 10% of the channel is deeper than 3.57 m, this result has important implications for how well a purely remote sensing approach can map the bathymetry of this stream, especially if deep pools are of interest.
For the Snake River, OPTID of image spectra led to OBRA R 2 values that stabilized beyond 1.8 m whereas field spectra led to OBRA R 2 values that increased up to the maximum depth in the data set: 1.9 m. Because the Snake was clearer and less deep than the Deschutes, the available data indicated that depths up to 1.8 m, and possibly deeper, were detectable by both the CASI sensor and the spectroradiometer but could not be used to constrain d max on the Snake. In this situation, knowing that depths up to at least a certain threshold can be resolved by a particular instrument is still useful information.
We also used irradiance profiles measured in the field to derive estimates of d max (λ) for the numerator and denominator wavelengths identified via OPTID of image and field spectra, indicated by point symbols in Figure 8b. A single value of d max could be obtained by taking the larger of d max (λ 1 ) and d max (λ 2 ) because even if one of the bands has saturated, such that further increases in d no longer lead to changes in R(λ), additional increases in d should still be detectable as long as the other band continues to change. Applying this approach to the Snake resulted in d max estimates that were much greater than the maximum depth observed in the field, implying that the CASI sensor could map the bathymetry of the entire reach without any gaps in pools or other deep areas.
An important advantage of using OPTID to estimate d max is that this empirical, image-based approach does not require specialized field equipment for making in situ optical measurements. Spectroradiometers are not widely available to river researchers and although we were able to use the same sensor to record both raft-based reflectance spectra and irradiance spectra at different depths within the water column, the latter involved a customized measurement crane and waterproof fiber-optic cable. Diffuse attenuation coefficients K d (λ) calculated from the resulting irradiance profiles effectively summarized the rate and spectral variation of attenuation of light by the water column in each river, which is the key physical process that enables depth retrieval from remotely sensed data. In order to calculate d max (λ) from K d (λ), however, we had to assume a value of ∆L N (λ)/L B (λ) as an index of detector sensitivity. How to specify this parameter a priori is not obvious; although various methods of estimating ∆L N (λ) have been proposed for oceanographic applications [49,50], these techniques are not readily applicable to fluvial environments. Identifying an appropriate value of ∆L N (λ)/L B (λ) thus remains ambiguous, making the direct, image-based OPTID approach an appealing alternative method for estimating d max .
The OPTID procedure provides a robust means of identifying the limitations associated with spectrally based bathymetric mapping and can thus help to promote more informed use of remotely sensed data. For example, an image-derived depth map based on OPTID output is presented in Figure 9a. The OBRA relation for the 3.57 m cutoff depth inferred to represent the sensor's d max was used to produce both the depth map and the calibration scatter plot shown in Figure 9c. This plot compares favorably with Figure 1a, which included all of the field data available for calibration. Truncating the data set increased the OBRA R 2 value and yielded a shallower minimum predicted depth of 0.25 m, as opposed to 0.73 m when all of the field data were used. The color scale for the depth map in Figure 9a extends from this minimum to the inferred d max of 3.57 m and thus provides a more detailed, realistic visualization of the plausible range of depth estimates than the bathymetric map in Figure 1b, which was produced without any kind of truncation or sampling. We also rearranged the X vs. d relation in Figure 9c to determine the X value corresponding to d max and then used this X value as an image threshold to mask out pixels likely to have depths in excess of d max . This threshold was exceeded for a very small proportion of the channel, less than 0.1%; the masked pixels are barely noticeable even in the zoom inset shown in Figure 9b. This representation indicates that the sensor's dynamic range for depth retrieval encompasses the majority of the channel. Moreover, truncating the range of input depths used for calibration improved depth retrieval within the shallow portions of the channel that are of primary interest for many practical, habitat-oriented applications. These results imply that even if spectrally based bathymetric maps are constrained by a maximum detectable depth, valuable information can still be derived from passive optical image data for those portions of the channel shallower than d max , which could be extensive in many cases. OPTID provides a rigorous approach to distinguishing-empirically, for a given image-between areas where depth estimates might be considered reliable and areas where depth estimates should be disregarded.

Extensions and Future Research
The sampling strategies presented herein were introduced in the context of hyperspectral images, but neither OPTID, SOBRA, nor the core OBRA algorithm necessarily require image data with high spectral resolution. As long as two or more bands are available, OBRA can be performed; even a basic red/green/blue (RGB) image provides three band ratios that can be compared to identify the optimal combination. Although hyperspectral data consisting of more bands open up the possibility of selecting narrow wavelength ranges highly responsive to changes in depth, Legleiter et al. [29] showed that resampling continuous field spectra to progressively broader bands only slightly decreased OBRA R 2 . Similarly, convolving the field spectra with the sensor response function of a multispectral satellite with only four visible and near-infrared bands lead to a marginal reduction in OBRA R 2 from 0.799 to 0.704. Moreover, two additional studies [24,31] produced realistic, reasonably accurate bathymetric maps from publicly available RGB aerial photography acquired on a regular cycle through the National Agricultural Imagery Program (NAIP). Previous research thus has shown that under appropriate environmental conditions (i.e., clear-flowing, relatively shallow rivers that are not obscured by vegetation), depth and reflectance are strongly related across broad portions of the visible spectrum.
Moreover, these results imply that bathymetry can be mapped accurately from standard image data sets acquired from aircraft, satellites, or, increasingly, unmanned aerial vehicles. The sampling strategies described in this study also can be applied to more widely available imagery of this kind. Because typical RGB images are recorded by using only 8 bits, limited radiometric resolution could constrain the dynamic range of image-derived depth estimates [36], making the ability to infer the maximum detectable depth via OPTID particularly valuable in an applied context.
This study introduced sampling strategies to improve spectrally based depth retrieval, but sampling is only one aspect of the overall process of calibrating a relationship between depth and reflectance. In addition to the data used for calibration, our focus herein, the method used to fit such a relation and the functional form of the equation or model used to link X to d are important components of the calibration process. We used ordinary least squares regression to fit second-order polynomials, but future studies should evaluate other statistical approaches and types of equations. For example, nonlinear and/or weighted least squares regressions could be considered along with exponential, power law, or logistic models. Improving upon all aspects of the calibration process-sampling, curve-fitting, and model selection-are important objectives to pursue if the potential benefits of remote sensing technology for river research and management are to be fully realized.

Conclusions
Although remote sensing techniques have demonstrated potential for characterizing river systems, obtaining reliable estimates of key parameters such as depth will require improved, more flexible methods and a better understanding of the uncertainties and limitations associated with a remote sensing approach. Motivated by this objective, we introduced a pair of depth retrieval algorithms that build upon the Optimal Band Ratio Analysis (OBRA) framework: OBRA of Progressively Truncated Input Depths (OPTID) and Stratified Optimal Band Ratio Analysis (SOBRA). The rationale for both procedures is to enhance bathymetric mapping performance by sampling field data to construct a calibration data set conducive to accurate depth retrieval across the full range of detectable depths in the river of interest. More specifically, OPTID creates a sequence of calibration data sets by iteratively discarding observations that exceed each one of a series of specified cutoff depths. SOBRA involves stratifying the data into a series of depth bins, each with the same number of observations, to avoid biases that might result from field methods that favor certain ranges of depths and thus lead to overor under-representation of some depth classes. Application of these methods to hyperspectral images and field spectra from the Deschutes and Snake rivers led to the following principal conclusions: 1.
OPTID can enhance bathymetric mapping by removing observations that exceed a sensor's maximum detectable depth from the calibration data set.

2.
Stratification can improve the accuracy of depth estimates relative to standard OBRA by ensuring that all ranges of depths are represented in equal proportions.

3.
Wavelength selection via OBRA is robust, can adapt to the distribution of depths included in the calibration data set, and tends to identify bands that are sensitive to variations in depth but do not saturate in deeper water.

4.
Identifying an appropriate cutoff depth to optimize depth retrieval yields insight regarding the range of depths for which reliable image-derived estimates can be obtained by a particular sensor deployed above a specific river.

5.
Retaining enough of the field data to enable sufficient variation in depth is also critical for establishing strong relationships between the image-derived quantity X and d; a trade-off must be reached between signal strength and saturation. 6.
OPTID can be used to estimate the maximum detectable depth directly from an image as an inflection point on a plot of OBRA R 2 vs. cutoff depth. This approach was consistent with d max (λ) values calculated from measured irradiance profiles on the Deschutes River but had the important advantages of not requiring specialized optical instrumentation or a priori assumptions of detector sensitivity. 7.
Although spectrally based bathymetric mapping is constrained by a maximum detectable depth, OPTID provides a robust means of distinguishing between shallow areas where depth estimates might be reliable and could provide useful habitat information and deeper areas where depth estimates should be disregarded.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/10/6/935/s1, Text S1: Text describing the effects of sampling strategy and sample size on Optimal Band Ratio Analysis (OBRA) of Progressively Truncated Input Depths (OPTID). The methods, results, and discussion associated with this analysis are all included in the Supplementary Information document, Figure S1: A figure summarizing the effects of sampling strategy and sample size on OPTID based on a hyperspectral image from the Deschutes River, Tables S1-S3: a set of tables summarizing the results of Stratified Optimal Band Ratio Analysis (SOBRA) based on a hyperspectral image from the Snake River and field spectra from both the Deschutes and Snake rivers. The tables included in the Supplementary Information document are analogous to Table 2