Assessing Stream Thermal Heterogeneity and Cold-Water Patches from UAV-Based Imagery: A Matter of Classification Methods and Metrics

Understanding stream thermal heterogeneity patterns is crucial to assess and manage river resilience in light of climate change. The dual acquisition of high-resolution thermal infrared (TIR) and red–green–blue-band (RGB) imagery from unmanned aerial vehicles (UAVs) allows for the identification and characterization of thermally differentiated patches (e.g., cold-water patches— CWPs). However, a lack of harmonized CWP classification metrics (patch size and temperature thresholds) makes comparisons across studies almost impossible. Based on an existing dual UAV imagery dataset (River Ovens, Australia), we present a semi-automatic supervised approach to classify key riverscape habitats and associated thermal properties at a pixel-scale accuracy, based on spectral properties. We selected five morphologically representative reaches to (i) illustrate and test our combined classification and thermal heterogeneity assessment method, (ii) assess the changes in CWP numbers and distribution with different metric definitions, and (iii) model how climatic predictions will affect thermal habitat suitability and connectivity of a cold-adapted fish species. Our method was successfully tested, showing mean thermal differences between shaded and sun-exposed fluvial mesohabitats of up to 0.62 ◦C. CWP metric definitions substantially changed the number and distance between identified CWPs, and they were strongly dependent on reach morphology. Warmer scenarios illustrated a decrease in suitable fish habitats, but reach-scale morphological complexity helped sustain such habitats. Overall, this study demonstrates the importance of method and metric definitions to enable spatio-temporal comparisons between stream thermal heterogeneity studies.


Introduction
The role of stream temperature as a key aspect of rivers' ecological functioning has received increasing attention in the last two decades [1][2][3]. Indeed, water temperature is a crucial player in all trophic levels and stages of biological organization in rivers. Water temperature governs the biochemical processes and conditioning of all freshwater life stages [4][5][6][7][8], highlighting the importance of understanding stream thermal heterogeneity. In addition, the predicted increase in water temperatures as a result of climate change [9] will likely grow the relevance of integrating the assessment and monitoring of stream temperature variability in future river management, particularly in identifying thermally resilient areas to changing climate [10][11][12][13][14].
Stream thermal heterogeneity is governed by a complex mosaic of interconnected habitats in rivers [15][16][17][18]. Though climatic, topographic, and hydrological controls drive stream temperatures at the air-water interface [1][2][3], multi-spatial scale factors such as in-than the mean reference, and their spatial dimensions can vary between 0.5 and 2.25 m 2 of minimum area [12,22,24,45,48,50,51]. Such variation in CWP threshold definitions combined with the complexity of UAV-based imagery data acquisition, processing, and analysis have reduced the comparability of thermal heterogeneity and CWP studies.
The identification of CWPs is a key first step towards understanding the role of potential stream thermal refugia and a valuable tool to support river conservation and sustainable management [52,53]. Predictions of CWP locations and distribution, for example, can help identify species competitive advantages or disadvantages over native species under climate change (e.g., [54]). The development of tools to predict the distribution of potential thermal refugia in a riverscape will be critical to understanding the resilience of biological communities in stream systems, to inform long-term measures that maintain resilient thermal landscapes in future climate-change scenarios [10,47,55], and to aid the development of effective climate adaptation tools [56].
With the aim to optimize approaches to assess thermal heterogeneity via UAV-based imaging and to help facilitate the comparability between stream CWPs research, the specific objectives of this study were: i.
To provide a standardized approach for the combined analysis of UAV-obtained RGB and TIR imagery, including a supervised classification approach of fluvial mesohabitats based on spectral properties (RGB), and a thermal heterogeneity assessment linked to the thermal properties (TIR) of the classified habitats. ii.
To quantify the changes in CWP numbers, distribution, and fluvial habitats with varying metric definitions, including thermal thresholds and patch size. iii.
To assess the relevance of thermal heterogeneity changes in future climatic scenarios, as well as the consequence of thermal habitat availability, using rainbow trout (Oncorhynchus mykiss) as a model species.

The Upper Ovens River
The Ovens River (Australia) is a 191-km-long north-east tributary of the Murray River with an annual discharge of 1776 GL.year −1 . The catchment area covers 0.7% of the Murray Darling basin and contributes 6% of the basin water [57]. The river upstream of its junction with the Buffalo River near Myrtleford is designated as the Upper Ovens River with an approximate length of 70 km [58]. The Upper Ovens has been highlighted as one of the most unaltered catchments within the Murray-Darling Basin, considered to be a high value system that provides suitable habitat conditions for many fish species, including the endangered Murray cod (Maccullochella peelii peelii), and the introduced cold stenothermic European brown trout (Salmo trutta) and North American rainbow trout (Oncorhynchus mykiss) [59]. Due to the highly connected groundwater and surface water resources in the Upper Ovens [19] and the ambient climatic conditions, the river provided an ideal study location for assessing stream thermal heterogeneity and the potential to provide CWPs for cold-water-adapted fish species during heat periods in the summer month through pockets of upwelling cooler groundwater [24].

Dataset and Representative Reaches
The dataset used for this study was obtained from the extensive UAV-based imagery survey presented by Casas-Mulet et al. [24], carried out in a 50 km long (linear length) reach of the Upper River at the end of summer 2017 through a total of 49 flights. Each flight consisted of simultaneously obtained RGB and TIR imagery across 580-1600 m river lengths. Single RGB images were acquired by means of Phantom 4 Pro drone (DJI, Shenzhen China) equipped with an optical red-green-blue band camera (Zenmuse X5 RGB, DJI, Shenzhen, China). Single TIR images were acquired by simultaneously flying a second drone (Inspire 1 Pro, DJI, Shenzhen, China) that was equipped with a thermal infrared camera (Zenmuse XT, DJI, Shenzhen, China). TIR images were calibrated using the ResearchIR MAX TIR imagery analysis software (FLIR systems Inc., Wilsonville OR, USA). The imagery set corresponding to each flight was merged into a single orthomosaic using the photogrammetry software Pix4D (Pix4D S.A., Prilly, Switzerland). To rectify basic overlay, we used field ground control points (GCPs) geo-located utilizing a Leica© VivaGS14 real-time kinetics (RTK) differential global positioning system (dGPS). The GCPs had been placed along the survey area, subject to the limitations of the riverscape. They were manually deployed before the flights so they could be visible in both TIR-and RGB-obtained imagery. The imagery's spatial resolution was 0.23 m per pixel in the TIR imagery and 0.035 m per pixel in the RGB imagery. For the TIR imagery, we used several in-situ water temperature point measurements and field-installed water loggers to correct and validate water temperature data, with a resulting correlation between logger value and corrected pixel data of R 2 = 0.7. We considered that the UAV-based TIR caption of surface water temperatures was representative of the measurements across the whole water depth due to the following reasons: (i) although our surveys were carried out during the lowest flow period, water was still running through the system so hydraulic mixing was ensured, and (ii) several spot temperature measurements in selected deep pools illustrated almost identical temperatures across the water column [60], which supported our assumption of no thermal stratification. We accounted for spatio-temporal variations in water temperature by acquiring the data over five days of stable and similar weather conditions in the summer, correcting daily variations (morning vs. afternoon) with the use of logger data, and considering altitude and latitude variations by using relative-to-reach average temperatures. Full details on data acquisition and processing can be found in the work of Casas-Mulet et al. [24].
For the analysis in this study, we selected five orthomosaics as representative reaches of contrasting morphology, riparian vegetation, shading input, and degree of disturbance. Each of the selected reaches represented the overall variability that repeats along the Upper Ovens riverscape. The selected reaches or survey sections (S1, S2, S3, S4, and S5) were considered representative for the system, as per the criteria defined by Braun et al. [16]. S1 and S2 (1210 and 810 m, respectively) represented the most unaltered and natural river sections, with broadly open and dynamic braided channels surrounded by natural riparian vegetation. S3 and S4 were more constrained than the previous sites because they were limited by geological factors; however, they were both still meandering channels (of 930 and 1600 m, respectively) mostly surrounded by agricultural land use and a thin strip of overhanging riparian vegetation. S5, of 580 m reach length, represented the most altered river section that was heavily channelized and straightened ( Figure 1). Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 21

Data Processing
To enable the combined thermal-physical habitat analysis of the Upper Ovens, we classified and combined both RGB and TIR imagery data into an accurate spatially matching dataset created using the World Geodetic System (WGS) 1984 system as a reference (Universal Transverse Mercator Zone 55S), through the geoinformatics software ArcMap 10.6.1 (Esri, Redlands, CA, USA). Such classification and dataset matching were essential for later data analysis, and they followed the four main steps illustrated in Figure 2 and described below: To enable the combined thermal-physical habitat analysis of the Upper Ovens, we classified and combined both RGB and TIR imagery data into an accurate spatially matching dataset created using the World Geodetic System (WGS) 1984 system as a reference (Universal Transverse Mercator Zone 55S), through the geoinformatics software ArcMap 10.6.1 (Esri, Redlands CA, USA). Such classification and dataset matching were essential for later data analysis, and they followed the four main steps illustrated in Figure 2 and described below: Figure 2. Functional block diagram illustrating the overview of both input and output data processes and analyses required in our methodological approach, going from drone-acquired RGB and thermal infrared (TIR) imagery to the resulting layer of intersection (LOI), followed by further data analysis.

Geo-Correction and Overlay
We geo-corrected each RGB mosaic to perfectly match the corresponding TIR imagery set and enable accurate overlay for further analysis. In addition to GCPs, we used highly visible points of contrasting temperatures (e.g., end tip of in-stream heated up deadwood) in both TIR and RGB mosaics as reference points for third-order polynomial correction. We used between 16 and 32 geo-correction points for each pair of mosaics, resulting in an overlapping RMS (root mean square) error of 0.28 m. We did not quantify the planimetric accuracy before and after adding both GCPs and highly visible points of contrasting temperatures in both TIR and RGB mosaics. However, we observed that the initial mosaicking using GCPs only was not enough to enable an accurate matching overlay between TIR and RGB despite the RTK dGPS accuracy. The most likely reason was their non-uniform and sparse distribution along the riverscape, which reduced planimetric accuracy [61,62]. GCPs had been manually deployed over a geometrically irregular area (longer than wider, which is characteristic of all river reaches), and their locations were limited to spots of riparian vegetation openings. Therefore, additional reference points of contrasting temperatures were essential to ensure an accurate matching overlay rather than being used as checkpoints. Functional block diagram illustrating the overview of both input and output data processes and analyses required in our methodological approach, going from drone-acquired RGB and thermal infrared (TIR) imagery to the resulting layer of intersection (LOI), followed by further data analysis.

Geo-Correction and Overlay
We geo-corrected each RGB mosaic to perfectly match the corresponding TIR imagery set and enable accurate overlay for further analysis. In addition to GCPs, we used highly visible points of contrasting temperatures (e.g., end tip of in-stream heated up deadwood) in both TIR and RGB mosaics as reference points for third-order polynomial correction. We used between 16 and 32 geo-correction points for each pair of mosaics, resulting in an overlapping RMS (root mean square) error of 0.28 m. We did not quantify the planimetric accuracy before and after adding both GCPs and highly visible points of contrasting temperatures in both TIR and RGB mosaics. However, we observed that the initial mosaicking using GCPs only was not enough to enable an accurate matching overlay between TIR and RGB despite the RTK dGPS accuracy. The most likely reason was their non-uniform and sparse distribution along the riverscape, which reduced planimetric accuracy [61,62]. GCPs had been manually deployed over a geometrically irregular area (longer than wider, which is characteristic of all river reaches), and their locations were limited to spots of riparian vegetation openings. Therefore, additional reference points of contrasting temperatures were essential to ensure an accurate matching overlay rather than being used as checkpoints.

Supervised Classification of RGB Imagery
For each RGB mosaic, we classified the key components of the riverscape, including fluvial mesohabitats (MHs), in-stream fluvial geomorphic features (GFs), and riparian cover (RC). MHs were classified into "run", "riffle", "pool", and "standing water" based on the appearance of the water surface structure to indicate flow velocity, depth visibility, and streambed substrate grain size [63][64][65]. We classified GFs into "gravel bar", "bedrock", and "dead wood". RC was classified into "riparian forest", "low vegetation", and "plantation". To account for the major effect that shading can have on stream temperatures [3], we distinguished between shaded and sun-exposed areas for each MH class and for the GFs classes "gravel bar" and "bedrock".
For the supervised classification of each RGB mosaic into the MH, GF, and RC classes, we performed a "maximum likelihood classification" (options tool intern setting at zero reject fraction and equal priori probability) in ArcMap. We first segmented each RGB mosaic into pixel sets of similar spectral properties. We then drew training polygons over the identified MHs, GFs, and RCs to sample the mean spectral properties for each class. (Figure 2). Each training polygon sample was converted into a reference signature file (data with mean RGB band values of each class) for each RGB mosaic. We used, on average, 142 training polygons per mosaic or survey section and 32 training polygons for each class. The resulting classified RGB mosaic was transformed into a polygon feature class (further referred to as RGB-classified layer; see Figure 2). To minimize automatically false classified pixel-segments, which were a result of inherent spectral variability [43], we eliminated polygons of true area in the field of ≤0.75 m 2 and manually assigned larger polygons that were wrongly classified to the correct class ( Figure 2).

Conversion of TIR to Homogeneous Thermal Patches
Given that thermal emissivity fluctuations from a water surface structure can influence temperatures assignments [66][67][68], the temperature values of the TIR imagery were rounded to the closest half-degree Celsius [49]. Each TIR mosaic was then transformed into a polygon feature class composed of multiple polygons of the same temperature (further referred to as thermal layer). Single polygons smaller ≤3 pixels (true area in the field of ≤0.16 m 2 ) were considered to be artefacts from the thermal camera sensor and were merged to the next larger polygon ( Figure 2).

Combining Classification and Thermal Layers into a Single Dataset
To enable the analysis of riverscape features and their temperature, the thermal layers were intersected with the corresponding RGB-classified layer, resulting in a final intersected layer (further referred to as layer of intersection-LOI). Since the TIR imagery pixel resolution was lower than the RGB imagery resolution, the intersection resulted in some ≤1 pixel (true area in the field of ≤0.053 m 2 ) polygons (or "silver polygons" [69]). These were eliminated by adding them to the next larger polygon of eachLOI ( Figure 2). The mean stream temperature of all MHs (true water areas) in each thermal layer was used as a reference for further analysis, and it was calculated by summing up the area values multiplied by the temperature values of each water polygon, then divided by the total water polygon area. The resulting mean stream temperature from the polygonised thermal layer was compared to the value obtained by averaging the water temperature pixel values from the original TIR raster layer. The results were the same, suggesting that the raster to polygon conversion and rounding did not affect the temperature precision.

Riverscape Classification Accuracy
To assess the reliability of the supervised classification method, we compared the classification outcomes to ground truth through a visual interpretation of the RGB imagery. We quantitatively assessed 30 sampling points per class, distributed randomly in each of the 5 RGB-classified layers. We used the Producer's accuracy indicator [70]  agreement of the supervised classification method with the visual interpretation, calculated as the percentage of agreement between predicted (classified) and observed (visual interpretation) data.
As an additional test of data reliability, we assessed the accuracy of land vs. water polygon classification in each LOI by looking at the discrepancy between the temperature values and MH classification of small (single pixels) polygons at the water-land edges. Basically, the MH polygons of single pixels with high-temperature values at the water-land edge are realistically gravel banks (not water) and should be classified as land. However, because of the non-perfect match between TIR and RGB mosaics, we found that some of these high-temperature pixels were classified as water. To assess to what extent this mismatch could produce thermal inaccuracies, we compared our reference mean water temperature calculation (see Section 2.3.4) to the same calculation that excluded these high-temperature 1-pixel size outliers for each LOI.

Stream Thermal Heterogeneity Assessment
For each LOI, we assessed the frequency (number) and area of occurring temperatures, relative to each mean reference temperature, in all MHs. Data were tested for normality using the Shapiro-Wilk test. Since the frequency data of MHs were not normally distributed, we used Kruskal-Wallis to test global significant differences. As a post-hoc test for individual pairwise comparisons of temperature distributions between sunny and shaded MH classes, we used the pairwise Wilcoxon rank sum tests with Bonferroni correction to adjust p-values for multiple testing.

Analysis of CWP Changes with Metric Definitions
We assessed changes in CWP numbers and distribution with temperature and area threshold definitions. For temperature, we considered a range between 1.5 and 3 • C (0.5 • C steps) cooler than the average stream temperature. We considered an ecological relevant minimum area size of 0.5 m 2 [22] as a starting point and increased it in steps to a maximum of 3 m 2 .
To assess the differences in CWP spatial distribution with changes in metric definitions, we computed the distances between the identified CWPs (start and end points being the start and end of the survey section) in each of three selected LOIs. These included one of each braided (S2), meandering (S4), and straight channel (S5). In order to ensure that distances were measured within the water body of the LOI, we used the attribute table tool in ArcMap to assign differentiated values between fluvial MHs (true water pixels), instream fluvial GFs, and RC classes and converted each LOI into a raster (connectivity raster; see Figure 2). Distances between CWPs were then analyzed with the "cost connectivity" tool in ArcMap using a total of 12 CWP metric definition combinations including areas of ≥1, ≥2, and ≥3 m 2 and relative temperature thresholds of ≤-1.5, ≤-2, ≤-2.5, and ≤−3 • C. To enable comparison between the selected LOIs or survey sections, the distances were calculated relative to each section's total water body length. CWP clusters (distances of <3 m between CWPs) were considered as a single CWP. Distances between CWPs grouped by metric definition and survey section were tested for normal distribution using the Shapiro-Wilk test. Global significant differences in CWP distribution between metric CWP definitions within each survey section were tested using the Kruskal-Wallis test.

Suitable Thermal Habitat Loss Assessment
We assessed how stream temperatures could change with an increased water temperature of 4 • C based on potentially predicted future climatic scenarios [71,72]. In order to assess how these changes in stream thermal heterogeneity would impact suitable thermal habitats for fish, we used rainbow trout (Oncorhynchus mykiss) as a target species. Rainbow trout is a temperature-sensitive fish species that is very relevant for fisheries globally and is particularly abundant in the mainstream Upper Ovens [59]. We considered 17 • C as a preference water Remote Sens. 2021, 13, 1379 9 of 19 temperature threshold for rainbow trout [73,74] and 22 • C as a maximum limit [61]. We used the same representative S2, S3, and S5 sections to illustrate the assessment.
All statistical analyses were carried out using the open-source statistical software R [75]. Significance was accepted at p ≤ 0.05.

Riverscape Classification Accuracy and Data Reliability
The supervised classification based on the spectral properties of water and their MHs, in-stream fluvial GFs, and RC resulted in a very high overall Producer's accuracy for all classes (Figure 3). MH classification resulted in slightly better results (mean = 85.4% and std = 10.22%) compared to the GF and RC classification (GF mean = 86.17% and std = 10.89%; RC mean = 80.46% and std = 17.12%). In considering each class separately, an overall higher variability in Producer's accuracy was found in shaded areas compared to sun-exposed MH (water) classes (Figure 3). Regarding land cover (GF and RC) classification, the lowest Producer's accuracy was related to the classes "deadwood" and "gravel bar shaded," whilst the highest Producer's accuracy was found for "bedrock" and "plantation" classes ( Figure 3). However, the "bedrock" class was only present in survey section S4. Our additional assessment of the effects of spatial mismatch that resulted in thermal inaccuracies suggested no significant differences between the reference mean water temperature calculation and the same calculation excluding high-temperature one-pixel size outliers for each LOI.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 21 globally and is particularly abundant in the mainstream Upper Ovens [59]. We considered 17 °C as a preference water temperature threshold for rainbow trout [73,74] and 22 °C as a maximum limit [61]. We used the same representative S2, S3, and S5 sections to illustrate the assessment. All statistical analyses were carried out using the open-source statistical software R [75]. Significance was accepted at p ≤ 0.05.

Riverscape Classification Accuracy and Data Reliability
The supervised classification based on the spectral properties of water and their MHs, in-stream fluvial GFs, and RC resulted in a very high overall Producer's accuracy for all classes (Figure 3). MH classification resulted in slightly better results (mean = 85.4% and std = 10.22%) compared to the GF and RC classification (GF mean = 86.17% and std = 10.89%; RC mean = 80.46% and std = 17.12%). In considering each class separately, an overall higher variability in Producer's accuracy was found in shaded areas compared to sun-exposed MH (water) classes ( Figure 3). Regarding land cover (GF and RC) classification, the lowest Producer's accuracy was related to the classes "deadwood" and "gravel bar shaded," whilst the highest Producer's accuracy was found for "bedrock" and "plantation" classes( Figure 3). However, the "bedrock" class was only present in survey section S4. Our additional assessment of the effects of spatial mismatch that resulted in thermal inaccuracies suggested no significant differences between the reference mean water temperature calculation and the same calculation excluding high-temperature one-pixel size outliers for each LOI.

Stream Thermal Heterogeneity Linked to Fluvial Mesohabitats (MHs)
The frequency and areal distribution of temperatures for each fluvial MH class followed non-normal distributions (Shapiro-Wilk test: p ≤ 0.05) (Figure 4), showing an overall heterogeneity without a particular trend or pattern by MH class. Stream temperature frequencies presented significant differences between MHs (Kruskal-Wallis test: p < 0.001). As expected, pairwise comparisons confirmed that sun-exposed MHs were significantly warmer than shaded MHs (pairwise Wilcox rank sum test: p < 0.001), with a maximum difference in temperature of 0.62 • C. Most MHs presented warmer temperatures than the average stream temperature, illustrating that cooler water areas were not dominant. The presence of MHs with relative temperatures of ≤−1.5 • C (lower CWP thermal threshold) was highly variable in most classes (between 1.49% and 6.47% in frequency and between 0.48% and 2.94% in area size), except for standing sun-exposed water types, where temperatures were always above this threshold.

CWP Types, Frequency, and Distribution Changes with Metrics
Changing CWP definition metrics within the established thermal and areal thresholds influenced the amount of identified CWPs, as well as the relative proportion of the fluvial MHs they covered ( Figure 5).
Changing CWP definition metrics within the established thermal and areal thresh-olds influenced the amount of identified CWPs, as well as the relative proportion of the fluvial MHs they covered ( Figure 5).
Following our expectation, the number of CWPs decreased as the thermal and areal thresholds were increased. Within a thermal threshold of ≤-1.5 °C, the number of CWPs decreased fivefold from ≥0.5 to ≥3 m²; however, the proportion of MH classes remained the same, following a similar pattern to the original river physical habitat composition (including both warm and cold patches). As thermal thresholds were increased, however, the proportion of MH classes changed from the original habitat composition and varied substantially between areal thresholds. The change in area share of MH composition was most sensitive with a shift from ≤-1.5 to ≤-2 °C. Such variable patterns, however, showed a consistent decrease in heterogeneity, with a loss of pool and riffle habitats, as area thresholds were increased ( Figure 5). The relative distribution of CWPs along each of the section types (braided, meandering, and straight) was not normally distributed (Shapiro-Wilk test p ≤ 0.05). Kruskal-Wallis tests revealed no global significant differences (p ≤ 0.14) in CWP distribution within each survey section. However, a visual assessment illustrated that changes in distribution were more sensitive to definition metrics in the straight channel (S5) than the braided channel (S2), which presented a smoother change in CWP distances based on metrics (Figure 6). Overall, the mean CWP distances increased with an increase of metric definition thresholds. However, the greatest increase in CWP distance occurred at different points of each of the representative reach-earlier (with lower temperature and areal thresholds) in the straight section and later (with higher temperature and areal thresholds) in the meandering and braided sections, suggesting that straight sections present faster changes and more complex sections present later shifts. Following our expectation, the number of CWPs decreased as the thermal and areal thresholds were increased. Within a thermal threshold of ≤−1.5 • C, the number of CWPs decreased fivefold from ≥0.5 to ≥3 m 2 ; however, the proportion of MH classes remained the same, following a similar pattern to the original river physical habitat composition (including both warm and cold patches). As thermal thresholds were increased, however, the proportion of MH classes changed from the original habitat composition and varied substantially between areal thresholds. The change in area share of MH composition was most sensitive with a shift from ≤−1.5 to ≤−2 • C. Such variable patterns, however, showed a consistent decrease in heterogeneity, with a loss of pool and riffle habitats, as area thresholds were increased ( Figure 5).
The relative distribution of CWPs along each of the section types (braided, meandering, and straight) was not normally distributed (Shapiro-Wilk test p ≤ 0.05). Kruskal-Wallis tests revealed no global significant differences (p ≤ 0.14) in CWP distribution within each survey section. However, a visual assessment illustrated that changes in distribution were more sensitive to definition metrics in the straight channel (S5) than the braided channel (S2), which presented a smoother change in CWP distances based on metrics ( Figure 6). Overall, the mean CWP distances increased with an increase of metric definition thresholds. However, the greatest increase in CWP distance occurred at different points of each of the representative reach-earlier (with lower temperature and areal thresholds) in the straight section and later (with higher temperature and areal thresholds) in the meandering and braided sections, suggesting that straight sections present faster changes and more complex sections present later shifts. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 21

Present and Future Availability of Suitable Thermal Habitat
A predicted increase in water temperatures illustrated a clear shift in thermally suitable habitats for rainbow trout (Figure 7). Habitats below the optimum thermal limit (≤22 °C) were still available to a large extent in the braided and meandering sections (99.9% and 98.4%, respectively) but reduced in the straight section (90.1%). Preferred thermal habitats (≤17 °C) were least available in the braided section (16.6%) and completely disappeared in both the meandering and straight sections. The latter was only left with unsuitable thermal habitats (Figure 7), likely resulting in the local extinction of the considered fish species.

Present and Future Availability of Suitable Thermal Habitat
A predicted increase in water temperatures illustrated a clear shift in thermally suitable habitats for rainbow trout (Figure 7). Habitats below the optimum thermal limit (≤22 • C) were still available to a large extent in the braided and meandering sections (99.9% and 98.4%, respectively) but reduced in the straight section (90.1%). Preferred thermal habitats (≤17 • C) were least available in the braided section (16.6%) and completely disappeared in both the meandering and straight sections. The latter was only left with unsuitable thermal habitats (Figure 7), likely resulting in the local extinction of the considered fish species.

Discussion
This study presents a reliable semi-automatic approach to extract both physical habitat and stream thermal temperature information from UAV-based RGB and TIR imagery and to enable their combined analysis to assess changes in thermal heterogeneity and CWPs. Our thermal heterogeneity assessment did not show a particular trend or pattern by MH class but illustrated significant differences between shaded and non-shaded types.

Discussion
This study presents a reliable semi-automatic approach to extract both physical habitat and stream thermal temperature information from UAV-based RGB and TIR imagery and to enable their combined analysis to assess changes in thermal heterogeneity and CWPs. Our thermal heterogeneity assessment did not show a particular trend or pattern by MH class but illustrated significant differences between shaded and non-shaded types. CWPs (using the lower relative temperature threshold of ≤−1.5 • C) were found in all except one MH; however, CWP numbers and the proportion of MHs they covered substantially changed with varying thermal and areal thresholds. The distribution assessment of CWPs revealed that this variability and metric sensitivity could be strongly linked to different types of reach morphology, with braided rivers being able to support bigger and cooler CWPs than straight channels. Such a reach morphology effect was also illustrated through the modelling of changes in a suitable thermal habitat with future predicted increases in water temperature. Channels with complex morphologies (braided and meandering) will provide more suitable thermal habitats for rainbow trout than simple straight channels in the future. These findings were based on the analysis of five characteristic reaches that reflected the full channel type variability along the riverscape. With each of them presenting a substantial reach length, they were considered representative of the ecohydrological processes occurring in the system. Overall, while acknowledging that any generalization needs to be made with care and that our dataset, although representative, was not exhaustive, we demonstrated that structurally homogenous river reaches are likely less resilient to increased temperatures than naturally braided ones in which CWPs are also sustained at increased ambient temperatures. Less complex or degraded channels can therefore bring potential local extinction consequences of cold-water-adapted species such as rainbow trout.

Riverscape Classification Accuracy and Data Reliability
Classification accuracy is crucial to assess the quality of thematic maps derived from remote sensing [70,76]. Our semi-automatic approach to riverscape classification was highly reliable to distinguish between riparian, in-stream, and water features, as well as to correctly classify them. A high Producer's accuracy of classification and a low difference in mean temperature calculation by excluding unrealistic high temperatures at the landwater interface attested to the high reliability and spatial accuracy of the generated dataset. Excluding temperature outliers in the LOI for analysis to increase the reliability of the dataset was more conservative, for example, than using a buffer zone around banks and tree canopy [45]. The latter approach ensures the exclusion of unintended terrestrial temperatures, but it was not suitable for our dual TIR and RGB imagery analysis because it could have excluded valuable water pixel data at the water-land edge representing CWPs provided by lateral or hyporheic discharge (e.g., upwelling zones). Our approach to separate water vs. land pixels made it possible to calculate mean stream temperature using all true water pixels, and it provided a comprehensive representation of the overall channel temperature. Our approach differs from the most commonly used ones in the literature, in which the calculation is often limited to selected points along a river thalweg [24,[45][46][47][48][49].
The fact that the data used in this study originated from two simultaneous but separated UAV flights (e.g., one for aerial RGB imagery and one for TIR imagery) posed a challenge for the exact georeferencing and overlapping of the two datasets. Still, our geocorrection approach achieved a precise matching between TIR and RGB with an RMS-error value similar to the smallest pixel size of the TIR imagery. Such a matching overlay enabled reliable water vs. land extraction and allowed for the analysis of stream water thermal heterogeneity. While our water vs. land extraction accuracy relied on RGB classification and was refined with the TIR-RGB discrepancy assessment (or LOI temperature outlier exclusion), such an approach was mostly possible due to the high contrast between land and water temperature in our dataset. In most cases, especially with low water-land contrast, the use of near infrared for water edge recognition would result in a higher accuracy of its delimitation. However, this system would still require an additional drone flight or camera mount, which was not available during our survey. Indeed, more advanced UAV setups that have two camera mounts or dual cameras combined with integrated, accurate spatial positioning technologies (e.g., RTK devices) may help the integration of multiple imagery data types in the future [77]. However, this recommendation can collide with regulations regarding UAV's maximum take-off mass and flight scenarios in the vicinity of urban areas or uninvolved people (e.g., EU drone regulations [78]), so our approach for processing and combing RGB-and TIR-derived imagery by two separate (lighter) UAV flights does not lose any relevance.

Stream Thermal Heterogeneity Links to Fluvial Mesohabitats (MHs)
The combined supervised classification method for RGB imagery with the accurate TIR imagery match resulted in a suitable approach to extract temperature values for each identified fluvial MH, helping overcome geo-correction and overlay challenges. Our classification of shaded vs. sun-exposed MHs confirmed that shading can have a major effect on cooling water temperatures [11,22]. The overall classification of CWPs revealed that cold water areas are not bound to specific MH classes and confirmed the findings by Casas-Mulet et al. [24], suggesting that CWP controls span over multiple spatial scales beyond the meso-scale.

CWP's Sensitivity to Metrics
Our analysis illustrated that the number, type, and distribution of CWPs change with varying definitions of temperature and areal thresholds. As expected, we found that as thermal and areal thresholds are increased, CWP numbers are reduced, their habitat heterogeneity decreases, and their distribution along the channel becomes sparser [10]. Changes in CWP numbers and distribution were found to strongly reflect different channel morphologies and landscape variables [12,24,50,79]. A single set of definition metrics provided significant differences in CWP distribution depending on the channel morphology; in straight channels, the number of CWPs decreased more rapidly with increased thermal and areal thresholds than in braided channels, reflecting the importance of complex morphology to sustain thermally differentiated areas. Regardless of the specific metric threshold, more complex river reaches can provide sufficient ranges of habitats that facilitate thermal heterogeneity and, in turn, promote CWPs; in contrast, straight channels require conservative thermal thresholds to enable the identification of differentiated thermal areas or CWPs due to their lower physical complexity [3,15,45,80,81]. The specific CWP definition metric thresholds used in this study fell within the range of definition metrics found in the literature [12,22,24,45,48,50,51,53], representing a broad variety of rivers in terms of size, morphology, and climatic context, as well as suggesting that our approach may enable potential comparisons. However, this study does not assume a direct transferability of results and does not attempt to provide a universal set of rules to identify and define CWPs. Rather, we consider our work the first development toward standard procedures and approaches for CWP assessments and comparisons.

Future Availability of Suitable Thermal Habitats for Fish
As expected, future climatic scenarios would reduce the spatial extent of thermally available habitats and reduce access to suitable refuge areas during climatic extremes. However, such reductions would be greater in less complex (i.e., structurally degraded) channels where thermally available habitats for fish would be substantially reduced as their distribution became sparser. Access to suitable habitats for rainbow trout was overall lower in simple straight channels than in complex braided reaches. Fish would hence need to swim further to reach suitable habitats in simple morphology reaches, suggesting a hindrance to fish population recruitment [82]. The structural deficits of river systems and their floodplains resulting from man-made channelization and river straightening are generally considered the main challenges in the conservation and restoration of riverine biodiversity (e.g., [13,83]). This is in line with our findings, highlighting the importance of in-stream river restoration measures to maintain thermally cooler water areas as potential thermal refugia in the future [84]. Our data, if computed at a large scale, might help to increase the number of remotely assessed water variables in the future. Future studies involving a combination of satellite and UAV-based TIR data can help predict large-scale thermal refugia occurrence, with a suite of potential applications to river ecologists and managers alike.

Conclusions
In conclusion, our approach to assessing combined UAV-based RGB-TIR imagery enabled an accurate integrated quantitative assessment of physical and thermal habitats. The identification of CWPs was sensitive to definition metrics but was also strongly linked to river morphological complexity, with complex channels able to support bigger and cooler CWPs in comparison to straight channels. Though this study does not attempt to provide a universal set of metrics to define CWPs, our findings help illustrate the importance of considering metric definitions and the fluvial geomorphology context to define CWP thresholds to (i) enable potential comparability between studies and (ii) develop standard procedures and approaches for CWP assessment in the future. The modelling of suitable thermal areas for cold-water fish estimated a decrease in suitable habitat for such species in the future, with warmer and less suitable waters compared to present conditions. The identified links with river morphological complexity and shading suggest that mitigation measures such as the restoration of channel complexity and establishment of bank vegetation can increase the resilience of fluvial systems to climatic change. The modelling approach presented here may generally be a powerful tool for predicting system resilience, prioritizing biological conservation, evaluating the effects of restoration measures, and supporting water resource management decisions.
In addition to its successful application in the Upper Ovens river, our approach can be easily transferred to other systems and climatic conditions, and results could be computed at a large scale using satellite data, for example. Our integrated approach should therefore be used to inform future climate adaptation measures and ensure the availability of suitable thermal habitats to sustain temperature-sensitive populations of freshwater biota.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.