Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An exploration of computational methods for classifying sediment patches within archived aerial photographs… Whitman, Peter 2019

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2019_september_whitman_peter.pdf [ 31.8MB ]
Metadata
JSON: 24-1.0380528.json
JSON-LD: 24-1.0380528-ld.json
RDF/XML (Pretty): 24-1.0380528-rdf.xml
RDF/JSON: 24-1.0380528-rdf.json
Turtle: 24-1.0380528-turtle.txt
N-Triples: 24-1.0380528-rdf-ntriples.txt
Original Record: 24-1.0380528-source.json
Full Text
24-1.0380528-fulltext.txt
Citation
24-1.0380528.ris

Full Text

AN EXPLORATION OF COMPUTATIONAL METHODS FOR CLASSIFYING SEDIMENT PATCHES WITHIN ARCHIVED AERIAL PHOTOGRAPHS OF GRAVEL-BED RIVERS by  Peter Whitman  B.A., Carthage College, 2017  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Geography)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  August 2019  © Peter Whitman, 2019   ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:  An Exploration of Computational Methods for Classifying Sediment Patches Within Archived Aerial Photographs of Gravel-Bed Rivers  submitted by Peter Whitman in partial fulfillment of the requirements for the degree of Master of Science in Geography  Examining Committee: Dr. Brian Klinkenberg, Geography Co-supervisor Dr. Michael Church, Geography Co-supervisor  Dr. Naomi Schwartz, Geography Additional Examiner    iii Abstract Bed material within gravel-bed rivers, which consists of gravel and other sediment coarser than 2mm, determines channel morphology and provides important benthic habitat. The impacts of flooding and anthropogenic activity on bed material are often examined to determine their effect on the morphology and ecology of gravel-bed rivers. To truly characterize this relationship, bed material must be discriminated from sand and sediment finer than 2mm at reach-scale over time. The current remote sensing routines that are used to synoptically characterize fluvial sediment did not emerge until the early 2000s, and as a result, reach-scale assessments of sand and gravel, extending beyond the 2000s, are absent for most gravel-bed rivers.  Fortunately, archived aerial photographs, can be used to analyze past landscapes. Traditionally, these analyses are carried out, manually, by photo interpreters, but multiple studies have shown that image processing techniques can be used to extract meaningful information from scanned aerial photographs. This study provides an exploration of semi-automated and automated image classification routines and their ability to replace manual interpretation for delineating patches of sand and gravel within scanned archived aerial photographs.   Results indicate that patches of sand and gravel within contemporary digital aerial imagery that has been degraded to mimic the characteristics of analog aerial photography can be consistently classified with overall accuracy above ~93% using automated object- and pixel-based classification routines. However, these same classifications only agree with manual interpretations of archived aerial photographs between ~45-70%. In contrast, the semi-automated routine provides measures of agreement that range almost entirely between ~75-85% when compared to the manual interpretations as well as the automated routines. Together, this demonstrates that a semi-automated routine should be used to classify scanned archived aerial photographs into sand and gravel.       iv Lay Summary Gravel-bed rivers are found throughout the world. The morphology of these rivers is largely determined by the erosion and deposition of gravel and sediment coarser than 2 mm. Retrospective characterization of gravel transport over time, especially as it relates to anthropogenic activity, requires that sand and sediment finer than 2 mm is discriminated from gravel. Archived aerial photographs of gravel-bed rivers act as tremendous historical record that can be manually classified into their sand and gravel components. However, this is a time consuming and potentially subjective task that requires highly skilled interpreters. The research herein compares multiple automated and semi-automated methods to determine if they can replace manual interpretation for classifying archived aerial photographs into their sand and gravel components. Results indicate that the semi-automated routine provides the best agreement with manual interpretation while producing a notably more detailed classification.    v Preface A portion of the ground-truthed data gathered in the field near Chilliwack, British Columbia in October 2017 was collected as part of a research project conducted by a graduate student—Chetan Chauhan-Sims—from Durham University in the United Kingdom and his supervisor—Dr. Patrice Carbonneau. However, the data have not been published yet by Chetan or Dr. Carbonneau.  Our research did not require ethics approval, nor has it been previously published in whole or in part. It was designed, carried out, and analyzed by the author of this manuscript and with advice from his supervisors.   vi Table of Contents Abstract ......................................................................................................................................... iii Lay Summary ................................................................................................................................ iv Preface ............................................................................................................................................ v Table of Contents .......................................................................................................................... vi List of Tables ................................................................................................................................. ix List of Figures ................................................................................................................................ x Acknowledgements ..................................................................................................................... xiii 1 Introduction ................................................................................................................................ 1 2 Literature Review ....................................................................................................................... 2 2.1 Gravel and Sand in Gravel-bed Rivers. ........................................................................... 2 2.2 Methods for Sampling Fluvial Sediment ......................................................................... 2 2.3 Aerial Photographs as a Historical Record of Past Landscape Dynamics ...................... 5 2.4 Summary .......................................................................................................................... 7 3 Research Objectives ................................................................................................................... 8 4 Methodology ................................................................................................................................ 9 4.1 Study Area ....................................................................................................................... 9 4.2 Ground-truth Data ......................................................................................................... 11 4.2.1 Data Collection ...................................................................................................... 11 4.2.2 Bootstrap Aggregated Generalized Linear Model ................................................. 12 4.3 Archived Aerial Photographs ........................................................................................ 15 4.3.1 Acquisition and Orthomosaic ................................................................................ 15 4.3.2 Resample ............................................................................................................... 17   vii 4.3.3 Segmentation ......................................................................................................... 19 4.4 Contemporary Aerial Images ........................................................................................ 19 4.4.1 Acquisition, Orthomosaic, and Segmentation ....................................................... 19 4.4.2 Degradation ........................................................................................................... 20 4.5 Classification ................................................................................................................. 23 4.5.1 Manual Interpretation ............................................................................................ 23 4.5.2 Object-based .......................................................................................................... 24 4.5.3 Pixel-based ............................................................................................................ 25 4.5.4 Semi-automated ..................................................................................................... 25 4.5.5 Agreement ............................................................................................................. 26 5 Results ........................................................................................................................................ 28 5.1 Ground-truth Model ....................................................................................................... 28 5.2 Sediment and Non-sediment Agreement ....................................................................... 30 5.3 Classification of Degraded Contemporary Orthomosaics ............................................. 31 5.4 Gravel and Sand Agreement .......................................................................................... 33 6 Discussion .................................................................................................................................. 38 6.1 Ground-truth Model ....................................................................................................... 38 6.2 Degradation of Contemporary Digital Aerial Imagery .................................................. 38 6.3 Classifications of Degraded Contemporary Digital Aerial Imagery ............................. 40 6.4 Agreement ..................................................................................................................... 41 6.4.1 Sediment and Non-sediment Classification ........................................................... 42 6.4.2 Gravel and Sand Classification ............................................................................. 44 6.5 Future Work ................................................................................................................... 48   viii 6.5.1 Ground-truth Model ............................................................................................... 48 6.5.2 Additional Variables for Classification ................................................................. 48 6.5.3 Degradation ........................................................................................................... 49 6.5.3.1 Granularity ......................................................................................................... 49 6.5.3.2 Generative Adversarial Neural Networks .......................................................... 49 6.5.3.3 Assessment ........................................................................................................ 50 7 Conclusion ................................................................................................................................. 51 References .................................................................................................................................... 53 Appendices ................................................................................................................................... 60 Appendix A ............................................................................................................................... 60 A.1 Ground-level Image Preprocessing ........................................................................... 60 A.2 Texture Statistics Generated from Grey-level Co-occurrence Matrix ....................... 62 A.3 Generalized Linear Model ......................................................................................... 64 Appendix B ................................................................................................................................ 69 B.1 Air Photo Preprocessing ............................................................................................ 69 B.2 Orthomosaic Workflow ............................................................................................. 72 B.3 Merge and Scale Parameters ..................................................................................... 74    ix List of Tables Table 4.1 The characteristics of the archived aerial photographs from each survey. .................... 16 Table 4.2 The ground sampling distance of each orthomosaic before and after they were resampled along with the nominal scale of each aerial survey ...................................................... 19 Table 4.3 The list of spectral, textural, and spatial attributes that are evaluated by the example-based feature extraction function in ENVI. ................................................................................... 26 Table 5.1 The output from the analysis of variance test. Each model is compared to the next most complex model to determine if there is a significant reduction in residual deviance when an additional independent variable is added to the model. ................................................................ 28 Table 5.2 The agreement between the semi-automated routine and the two photo interpreters when classifying an example archived aerial photograph from each survey into sediment and non-sediment. ................................................................................................................................ 31 Table 5.3 The optimal merge and scale parameters for each aerial survey. .................................. 31 Table 5.4 Agreement between the two photo interpreters, semi-automated routine, and automated pixel- and object-based routines when classifying an example archived aerial photograph from each survey into sand and gravel. .................................................................................................. 35    x List of Figures Figure 4.1 An outline of the methodology with section 4.2 in green, 4.3 in blue, 4.4 in orange, and 4.5 in purple. ............................................................................................................................. 9 Figure 4.2 A map of the lower course of the Fraser River (from Church, 2017). ......................... 10 Figure 4.3 An illustration of the method used to collect ground-level image data (from Graham et al., 2005). ....................................................................................................................................... 11 Figure 4.4 An example of a training image (left) and the corresponding binary representation of its sand and gravel fractions produced using Adobe Photoshop (right). ....................................... 12 Figure 4.5 A skewness-kurtosis plot of the empirical distribution along with common theoretical distributions produced using the fitdistrplus package in R (Delignette-Muller, Dutang, Pouillot, Denis, & Siberchicot, 2019). ......................................................................................................... 13 Figure 4.6 A conceptual framework of bootstrap aggregation. ..................................................... 15 Figure 4.7 The workflow used to georeference each orthomosaic derived from archived aerial photographs. .................................................................................................................................. 17 Figure 4.8 An example of a sample collected perpendicular to an edge (top) the resulting edge spread function (bottom-left) and the line spread function (bottom-right) produced by taking the derivative of the edge spread function. .......................................................................................... 18 Figure 4.9 The red line represents the ~13km flight used to collect 12 digital aerial images in March 2018. The blue circles denote the approximate principle point of each image while the yellow boundary represents the areal extent of the orthomosaic derived from the images. .......... 20 Figure 4.10 The spectral sensitivity of Kodak Double-X 2405 aerographic film (from Eastman Kodak Company, 2005), with an orange dashed line representing the 525 nm filter used in the 1984, 1999, and 2008 aerial surveys. ............................................................................................ 21   xi Figure 4.11 The spectral sensitivity of Vexcel Ultracam Eagle large format digital camera (from Vexcel Imaging, 2015), with an orange dashed line representing the 525 nm filter used in the 1984, 1999, and 2008 aerial surveys. ............................................................................................ 22 Figure 4.12 A schematic of the histogram matching transformation. ........................................... 23 Figure 4.13 The scheme used to categorize the manually interpreted classes. ............................. 24 Figure 5.1 The predictions of percent gravel produced by the chosen model plotted against observed percent gravel. For comparison, a 1:1 line is provided (dotted orange) along with a locally weighted smooth line (solid blue). .................................................................................... 29 Figure 5.2 A comparison of mean and median as aggregators and the number of samples required to stabilize root mean square error. ................................................................................................ 30 Figure 5.3 The overall accuracy of 100 bootstrapped objected-based classifications for each degraded contemporary orthomosaic. ........................................................................................... 32 Figure 5.4 The overall accuracy of 100 bootstrapped pixel-based classifications for each degraded contemporary orthomosaic. ........................................................................................... 33 Figure 5.5 The agreement between 100 bootstrapped object-based (top) and pixel-based (bottom) classifications when applied to the sediment regions of an example archived aerial photograph from each survey ........................................................................................................................... 37 Figure 6.1 Points with ground-truthed estimates of percent gravel along with the contemporary orthomosaic after its degradation to match the aerial photographs from 2008. Each point represents a “site”, which is the mean center of four samples of percent gravel taken within a visually homogenous patch of sediment. The red label at each site location provides the average percent gravel value for the surrounding sediment patch. ............................................................. 40   xii Figure 6.2 The outline of the four contact prints used to assess agreement between various classification schemes. .................................................................................................................. 41 Figure 6.3 Examples of sediment/non-sediment classification: original (top), manual interpretation (bottom left), semi-automated routine (bottom right). Non-sediment regions are in black. ............................................................................................................................................. 43 Figure 6.4 Examples of the two interpreter’s classification of the 2008 contact print into its sand and gravel components. Non-sediment regions are in black, sand is denoted by the colored regions, and gravel makes up the remaining areas of each image. ................................................ 44 Figure 6.5 An example of each automated routine used in this study: object-based SVM classification (top left), object-based RTC classification (top right), pixel-based SVM classification (bottom left), pixel-based RTC classification (bottom right). ................................. 46 Figure 6.6 An example of the semi-automated classification of the 2008 contact print into its sand and gravel components. Non-sediment regions are in black, sand is denoted in orange, and gravel makes up the remaining areas of the image. ....................................................................... 47 Figure 6.7 The primary ordering of image analysis elements that were used by each classification paradigm in this study (adapted from Estes, Hajic, & Tinney, 1983). .......................................... 49     xiii Acknowledgements First and foremost, I would like to thank my co-supervisors, Dr. Brian Klinkenberg and Dr. Michael Church for their patience and guidance in all aspects of my graduate program. They were able to draw on a tremendous amount of past experience to perfectly balance autonomy with prompt assistance when necessary. Funding was provided by operating grants to Dr. Church from the Natural Sciences and Engineering Research Council (NSERC) of Canada. The author was also supported by a Graduate Support Initiative Award from the UBC Faculty of Arts. Thank you to my lab mates Emily Acheson, Mielle Michaux, and Michael Jerowsky who were present during the glamorous and unglamorous moments of graduate school. Each of them exhibit qualities that I admire immensely and try to emulate in and outside of academia. Additional thanks to my father, mother, brother, and grandfather who inspired me to pursue a graduate degree and provided continuous advice, support, and love throughout my master’s program.      1 1 Introduction  Gravel-bed rivers shape the upper reaches of fluvial systems throughout the world. The bed material within these rivers, which consists of gravel and other sediment coarser than 2mm, dictates channel morphology and provides important benthic habitat. There has been ongoing work to characterize the relationship between bed material and flooding as well as anthropogenic activity to determine how the morphology and ecology of gravel-bed rivers is impacted by exterior forces. To establish this relationship, bed material must be examined over decadal scales, and in doing so, differentiated from sand and sediment finer than 2mm (Ham, 2005). This is a difficult task because rivers are dynamic and continuous systems that vary temporally and spatially throughout the length and scale of an entire reach (Marcus & Fonstad, 2008). Furthermore, remote sensing routines, which are best adapted to accommodate reach-scale variability in fluvial sediment, did not emerge until the early 2000s. Before that, sediment was monitored using ground-level techniques that are limited in their spatial and temporal coverage. Consequently, reach-scale assessments of sand and gravel before the 2000s are absent for most gravel-bed rivers.  Fortunately, archived aerial photographs, collected as early as the 1920s in many parts of the world, act as a robust historical record of past landscapes. While aerial photographs are traditionally manually classified by photo interpreters, multiple studies have used image processing techniques on scanned aerial photographs to extract meaningful information regarding topography and ecological dynamics (Cots-Folch, Aitkenhead, & Martinez-Casasnovas, 2007; Lane et al., 2010; Morgan & Gergel, 2010; Okeke & Karnieli, 2006a, 2006b). This study aims to explore prominent semi-automated and automated image classification routines and their ability to replace manual interpretation for delineating patches of sand and gravel within scanned archived aerial photographs of gravel-bed rivers. Archived aerial photographs are usually devoid of reference data that can be used to validate classification efforts. Therefore, the automated routines, in particular, will be explored to determine if a classification can be developed and validated using ground-truth data paired with contemporary digital aerial images after they have been degraded to mimic the characteristics of archived aerial photographs.   2 2 Literature Review 2.1 Gravel and Sand in Gravel-bed Rivers. The form and process of a river is largely dictated by the erosion, transportation, and deposition of sediment (Church, 2006). The diameter of individual clasts within sediment can range nominally from clay, silt, and sand to gravel, cobbles, and boulders (Wentworth, 1922). The sediment dimensions associated with these classes influence channel morphology and hydraulics, as well as fluvial ecosystem function (Rice & Church, 2010).  Fluvial sediment within a defined reach of a river is often categorized into bed material and wash material. Bed material determines the morphology of the reach and makes up the channel bed and lower banks. In contrast, wash material passes through a reach once entrained and is deposited only in quiet backwaters and overbank during floods. As the name suggests, in gravel-bed rivers, bed material consists of gravel and other sediment coarser than ~2 mm and wash material consists of sand and other sediment finer than ~2 mm (Church, 2006). From here on, bed and wash material will be loosely equated to gravel and sand, respectively.    2.2 Methods for Sampling Fluvial Sediment  Most transported bed material in gravel-bed rivers is stored in large-scale sediment bars (Neill, 1987). When exposed during low flow, these bars allow researchers to measure grain size distributions, which provides information about sediment load that can inform river management practices and sediment budget calculations. However, due to the dynamic and continuously varying nature of fluvial systems, even visually homogeneous patches of sediment along mid-channel bars can exhibit enough variation in grain size to necessitate robust sampling methods (Buffington & Montgomery, 1999; P. Nelson, Bellugi, & Dietrich, 2014). The traditional approach to calculating surface grain size distributions is the grid-by-number technique, which involves blindly selecting and then measuring grains from the bed surface in a grid pattern (Wolman, 1954). An investigation into the errors associated with this technique found that percentile precision improved with an increasing number of samples (Rice & Church, 1996). While smaller samples can be used to determine the median grain size, to produce a representative grain size distribution with a reasonably strict confidence interval Rice & Church (1996) recommend that samples consist of 400 stones. This sampling approach, along   3 with other established manual, field-based techniques, is accurately described as onerous. They require extended effort in the field, disturb the sediment substrate, and lead to samples with limited spatial coverage.  Ground-level image processing techniques for grain size measurements first emerged in the 1980s. These techniques are carried out by taking images of sediment either in the field or during a flume experiment and then processing those images in such a way that yields grain size measurements. Two main categories of techniques exist: direct or object-based approaches and indirect or statistical approaches (Black, Carbonneau, Church, & Warburton, 2014). Object-based approaches measure individual particles within an image while statistical approaches relate grain size to measures of spatial variation in pixel intensity values. While they don’t completely solve the problem of limited spatial coverage, ground-level image processing techniques quickly gained popularity because they are less invasive, reduce the time required in the field, and provide the potential for automation.   One of the first examples of ground-level image processing is an object-based approach published by Ibbeken and Schleyer (1986) that involves measuring the axial lengths and orientation of ellipsoids manually fit to the outline of each particle in an image. More recently, object-based approaches have been automated. Graham, Reid, and Rice (2005) explored several types of image segmentation routines to isolate particles within an image so that ellipsoids could be fit automatically to those particles for subsequent measurement of axial length and orientation. This work was incorporated into a software called BASEGRAIN, which is used widely within hydrological and morphological studies (Detert & Weitbrecht, 2012).  Studies investigating statistical-based approaches started to emerge in the early 2000s. These studies show that measures of texture, computed for an image using spatial autocorrelation or a grey-level co-occurrence matrix (GLCM), vary with grain size (Carbonneau, Bergeron, & Lane, 2005; Rubin, 2004). Spatial autocorrelation is a measure of the correlation between the intensity values of each pixel in a neighborhood of an image (Cliff & Ord, 1973). Rubin (2004) presented a methodology in which measures of spatial autocorrelation are calibrated to grain size distributions produced by traditional ground-based measurements. Similar to spatial autocorrelation, GLCMs measure the distribution of co-occurring pixel values within a neighborhood of an image (Haralick, 1979; Haralick, Shanmugan, & Dinstein, 1973). They can   4 be used to calculate multiple statistics such as entropy, contrast, homogeneity, correlation, energy, and dissimilarity (Clausi, 2002; Haralick, 1979; Haralick et al., 1973). Carbonneau et al. (2005) used entropy derived from a GLCM to segment ground-level digital images into their constituent sand and gravel fractions.  Despite the advances provided by the outlined image processing methods, all ground-level techniques require that samples be taken at the bar surface, which results in samples with limited spatial coverage. Remote sensing routines involving digital aerial imagery are an alternative to ground-based techniques that allow sediment to be sampled continuously at reach-scale. Statistical methods are primarily used to carry out these assessments because images collected by aerial platforms are of lower spatial resolution than images captured at ground level. Pixel values are a function of the average spectral response for the area that the pixel covers. In the case of images taken at ground-level, pixel size is much smaller than grain size, which makes it possible to delineate individual grains using object-based methods. However, most reach-scale assessments involve aerial images captured by airplanes that exhibit spatial resolutions of 3 cm to 18 cm (Black et al., 2014; Carbonneau, Lane, & Bergeron, 2004; Chandler, Rice, & Church, 2004; Verdú, Batalla, & Martínez-Casasnovas, 2005). The pixels within these images are often large enough to encompass multiple grains, making it difficult to delineate individual particles.  In a preliminary study, Chandler et al. (2004) attempted to classify aerial images at multiple scales into one of five discrete classes: sand, sandy gravel, pebble, clean gravel, and cobble. The authors used a supervised classification that was trained and validated with georeferenced field-based measurements of grain size and found that the overall accuracy of their predictions improved once texture, computed using a standard deviation filter, was included in their classification scheme. Surprisingly, the authors also discovered that scale did not appear to have a significant effect on their results.  Carbonneau et al. (2004) and Verdú et al. (2005) were able to map continuous estimates of grain size at reach-scale by developing predicative models based on the relationship between georeferenced ground-level grain size measurements and image texture. Similar to many ground-level image processing methods, both studies quantified texture using spatial autocorrelation and statistics derived from a GLCM. The predictive relationships were validated using independent data collected through ground-level manual sampling. Both Carbonneau et al. (2004) and Verdú   5 et al. (2005) tested their methodologies using images collected at multiple scales. They discovered, contrary to Chandler et al. (2004), that the accuracy of their predictions increased when using images with higher spatial resolution.   Black et al. (2014) blended the approaches outlined in previous studies to analyze aerial imagery with 3 cm spatial resolution. Similar to Chandler et al. (2004), the authors delineated sand patches from coarser patches of sediment by applying a binary segmentation routine to texture images produced using a standard deviation filter. The resulting overall accuracy was ~91%, which is much higher than the accuracy achieved by Chandler et al. (2004). This is largely due to the number of classes delineated within each study. Chandler et al. (2004) attempted to produce five classes and found that their classification performed better for cobble (~81%) than for the classes representing smaller clasts such as sand (~39%) and sandy gravel (~43%). One assumes that, if Chandler et al. (2004) had reduced the number of classes, their overall accuracy would have increased. Inspired by Carbonneau et al. (2004) and Verdú et al. (2005), once the sand and gravel fractions were delineated within their images, Black et al. (2014) mapped grain size within the coarse fraction by developing predictive models predicated on multiple GLCM statistics.  The studies outlining statistical methods for digital aerial images as well as ground-level images show that there is a strong relationship between grain size and measures of image texture regardless of image scale. Black et. al (2014) did highlight that this relationship is strongest when pixel size is smaller or larger than grain size and breaks down when spatial resolution and grain size are equal. Black et. al (2014) also noted that multispectral images are not needed to extract meaningful information about grain size because all wavelengths accounted for a similarly high level of variance in spectral intensity within sediment regions. It is therefore conceivable that these fundamental aspects of analyzing fluvial sediment using digital images could be exploited to analyze archived aerial photographs once they are scanned.  2.3 Aerial Photographs as a Historical Record of Past Landscape Dynamics Aerial photographs were captured using film cameras. Instead of a sensor that records spectral information as a matrix of intensity values, this information is recorded on negative film using a photosensitive emulsion that is usually comprised of silver halide crystals suspended in   6 gelatin. These crystals darken when exposed to light, which results in a negative that represents the reverse of what the human eye perceives. After the negative film has been exposed, it is developed to produce a positive enlargement or contact print. During this process, the negative film goes through a reversal so that the resulting positive displays light and dark areas as expected. Aerial photographs tend to be either panchromatic, true color (RGB), near-infrared (NIR), or color-infrared (CIR) with spatial resolution depending largely on flying height and the photo development process (Morgan, Gergel, & Coops, 2010). In Canada and many other parts of the world, systematic acquisition of aerial photographs began as early as the 1920s and 30s. Even though they have largely been replaced by digital aerial images for most purposes, archived aerial photographs remain a tremendous source of information for studies focused on past landscapes (Morgan et al., 2010). Examples of these studies include efforts by ecologists to quantify changes in land cover and landscape heterogeneity (Cots-Folch et al., 2007; Morgan & Gergel, 2010; Okeke & Karnieli, 2006a, 2006b), and by fluvial geomorphologists to track changes in surface elevation (Lane et al., 2010). For many years, analysis of aerial photographs was limited to manual interpretation. Analog photographs are now commonly scanned so that they can be analyzed as digital images. This is often considered more objective, consistent, and less time consuming than traditional photo interpretation (Brasington, Langham, & Rumsby, 2003; Morgan et al., 2010). Despite these advancements, analyses involving archived aerial photographs are difficult to validate because georeferenced data are usually nonexistent. Building on the work of Kelly (1960) surrounding the quantification of the photographic development process, Geigel & Musgrave (1997) provide an outline of the steps needed to simulate the attributes of analog photographs in digital images. Specifically, they identified exposure, spectral sensitivity, resolving power, and granularity as the four primary attributes that must be emulated to degrade digital images so that they become facsimiles of a target analog photograph. So far, the steps outlined by Geigel & Musgrave (1997) have only been used to degrade hand-held digital images, but they could potentially be used to alter contemporary aerial images so that they become facsimiles of archived aerial photographs. If these degraded aerial images are paired with georeferenced data, a classification scheme could be developed, validated, and then used to analyze archived aerial photographs with more confidence.   7 2.4 Summary There are remote sensing routines involving contemporary digital aerial images that can be used to map fluvial sediment at reach-scale. These methodologies indicate that there is a strong relationship between image texture and grain size that could be exploited in scanned archived aerial photographs to delineate patches of sand and gravel. However, as it stands, any analyses conducted with aerial photographs are going to be devoid of validation data. Consequently, there is additional motivation to discern if a classification routine can be developed and validated using degraded contemporary digital aerial images and applied to archived aerial photographs.     8 3 Research Objectives The main aim of this study is to investigate the ability for image processing techniques to produce a method capable of delineating patches of wash material (sand; less than or equal to 2mm) and bed material (gravel; greater than 2mm) within archived aerial photography.  Drawing from this aim, the following research objectives have been identified: 1. Develop a flexible model to determine the percent sand and percent gravel within ground-level digital images of sediment bar surfaces 2. Simulate the attributes of archived aerial photographs in contemporary aerial images, so that the two sets of remotely sensed products can be analyzed as facsimiles of one another.   3. Compare pixel-based and object-based classification schemes in order to determine the most accurate technique for mapping sand and gravel patches in the degraded contemporary aerial images. 4. Formulate a reproducible method capable of mapping surface sediment patches for archived aerial photography using the optimal classification technique developed with the degraded contemporary images.   9 4 Methodology This section will begin with a description of the study area (4.1) and then follow the general workflow provided in Figure 4.1 with subsections reviewing methods associated with preparing and analyzing field data (4.2), archived aerial photographs (4.3), contemporary digital aerial images (4.4), and image classifications (4.5).   Figure 4.1 An outline of the methodology with section 4.2 in green, 4.3 in blue, 4.4 in orange, and 4.5 in purple.  4.1 Study Area The Fraser River drainage basin covers an approximate area of 232,000 km2, and stretches nearly 1400 km from its headwaters near the British Columbia and Alberta border to its delta on the Pacific Ocean. At Hope townsite, ~1000 km downstream from its headwaters, the river transforms into its lower course and flows through a broad alluvial plain where the gradient declines, sediment is deposited, and the channel becomes unconfined. This study is concerned with the 55 km gravel reach that extends from Laidlaw to Sumas Mountain (Figure 4.2). As the   10 name suggests, the gradient within this reach has declined enough that gravel and coarser sediments are deposited and stored as bed material in a complex system of mid-channel bars and islands. Sand and finer sediments, on the other hand, mostly remain entrained as they travel downstream to the sand-bed reach (Venditti & Church, 2014), but they can be transiently deposited on bar tops and in back channels within the gravel-bed reach.   Figure 4.2 A map of the lower course of the Fraser River (from Church, 2017).  Annually, the flow of water through the gravel reach exhibits a unimodal pattern. The highest monthly average discharge measured at Hope townsite occurs in June as a result of the spring freshet (Environment Canada, 2016). Discharge then returns to a relatively stable low flow that extends from October to April. During peak discharge, the water level rises an average of ~4 meters, which puts the surrounding area at risk of flooding and inundates the mid-channel bars and islands within the reach. This annual process can contribute to coarse sediment transport and form deposits of fine sediments eroded from channel banks along the mid-channel bars and islands. Collectively, the magnitude of these changes depends on the severity and frequency of flooding, but can take several years to detect (Ham, 2005).       11 4.2 Ground-truth Data 4.2.1 Data Collection  Fieldwork was conducted during a six-day campaign in October 2017 and a one-day campaign in March 2018 to collect georeferenced images of surface sediment at 191 sites along the channel bars of the gravel reach. These sites were located within visually distinct patches of sediment that ranged continuously in composition from homogenous sand to homogenous gravel. Four RGB images of the bar surface were taken at random locations surrounding each site using a hand-held digital camera (Figure 4.3). The locations of these images were then recorded using a Garmin eTrex 20x GPS unit that provided an estimated accuracy of ± 5m (Garmin, 2011). Accordingly, all patches chosen as survey sites had minor axes greater than or equal to ~10m so that the recorded image locations would not be confused with adjacent sites.   Figure 4.3 An illustration of the method used to collect ground-level image data (from Graham et al., 2005).  Each image captured in the field contained a 0.75m x 1.0m wooden frame so that they could be scaled. In an initial step, the images were cropped by the interior dimensions of that frame to ensure that its presence would not confound any subsequent analysis. The MATLAB   12 code for this automated process and a visual representation are provided in Appendix A.1. After they were cropped, the images were visually sorted into one of three categories: homogenous sand, homogenous gravel, and mixed. With ~200 images containing a mixture of sand and gravel, a predictive model was developed to avoid the time consuming and potentially subjective task of visually estimating the percent cover.  4.2.2 Bootstrap Aggregated Generalized Linear Model A generalized linear model was developed to estimate the percent sand and percent gravel within each mixed image by exploiting the established relationship between image texture statistics and grain size. This model was trained using 30 images that were split into quadrants to artificially boost sample size. Together, the resulting 120 images provided the full spectrum of variation in gravel and sand cover. The response variable within the model is the percent cover of gravel in an image and is represented by 0	 < 	$	 < 	1. This variable was provided for each training image by manually delineating the gravel fraction using Adobe Photoshop. A binary image was produced, which allowed easy estimation of sand and gravel cover as a percentage of total image area (Figure 4.4).    Figure 4.4 An example of a training image (left) and the corresponding binary representation of its sand and gravel fractions produced using Adobe Photoshop (right).  Following Cullen & Frey (1999), an appropriate distribution to model the response variable was selected by plotting the skewness and kurtosis characteristics of the empirical distribution along with common theoretical distributions (Figure 4.5). The characteristics of 1000   13 bootstrapped samples were also plotted to account for uncertainty. The resulting plot shows that the response variable is best modeled using the beta distribution.   Figure 4.5 A skewness-kurtosis plot of the empirical distribution along with common theoretical distributions produced using the fitdistrplus package in R (Delignette-Muller, Dutang, Pouillot, Denis, & Siberchicot, 2019).  Six independent variables were generated for greyscale versions of each training image: brightness, contrast, correlation, energy, entropy, and homogeneity. Aside from brightness, which represents the average intensity value within each grey-scale image, all of the independent variables are measures of image texture. These statistics were generated from a grey-level co-occurrence matrix (GLCM) (Appendix A.2), which calculates how often pairs of pixels occur in an image (Haralick et al., 1973). Following the recommendations of Clausi (2002), the computational expense of generating the GLCMs was reduced by scaling all of the images from 256 grey-levels (8bit or 28) to 8 grey-levels (3bit or 23).    14 Candidate models were constructed by generating all possible combinations of independent variables up to sets of three. These models were also given a random effect that indicates the day that the image was taken during the field campaigns to account for variance in lighting conditions. Including a null model, this produced 42 candidate models that can be expressed using the following equation:  $&' = (β+	+	-+') + (β/	+	-/')x/&' +	. . .		(23 +	-3')43&' + ϵ&'  where $&' denotes the percent gravel in the i-th image of day j, 46&' and 7&' are the respective independent variables and error term, and -6' the random slope coefficients for group j. To select the best performing model, candidate models with the lowest Akaike information criterion (AIC) value at each unique level of degrees of freedom were compared using an analysis of variance (ANOVA) test.  Once the optimal model was chosen, bootstrap aggregating (bagging) was employed to improve accuracy, reduce variance, and avoid overfitting (Figure 4.6). Bagging is predicated on the idea that the predictions of multiple weak models aggregated together outperform the predictions of a single strong model (Breiman, 1996). The type of aggregator—mean or median—and the number of bootstrapped samples needed to outperform a single strong model were examined by computing the RMSE after each new sample had been aggregated.     15  Figure 4.6 A conceptual framework of bootstrap aggregation.   After estimates of percent gravel had been computed for each mixed image, they were combined with the separated samples consisting of homogenous sand or homogenous gravel to form a single data set of percent gravel estimates at each sample location. The geographic mean center of the constituent samples within each site was computed using ArcGIS 10.5.1 (Environmental Systems Research Institute, 2016). The average estimate of percent gravel was appended to that location and classified into sand (< 50% gravel) and gravel (≥ 50% gravel).  These georeferenced locations of sand and gravel were then used in subsequent methodological steps as ground-truthed calibration and validation data.  4.3 Archived Aerial Photographs 4.3.1 Acquisition and Orthomosaic A set of archived aerial photographs from aerial surveys that occurred in 1951, 1984, 1999, and 2008 were acquired from the Geographic Information Center (GIC) at the University of British Columbia (UBC) and the National Earth Observation Data Framework Catalogue (NEODF-Cat), which is maintained by Natural Resources Canada in Ottawa. These particular sets of aerial photographs were chosen because they correspond to hydrographic surveys that were carried out in those years. This makes them convenient for application to the ongoing analysis and improvement of sediment budget estimates for the Fraser River. All of the photographs were captured on panchromatic film and acquired as 10in x 10in contact prints   16 (Table 4.1). These photographs were scanned at 2400 dpi to create TIFF images with 256 grey-levels (8bit) using an EPSON Expression 120000XL scanner.   Table 4.1 The characteristics of the archived aerial photographs from each survey.  1951 1984 1999 2008 Date 1951-OCT 1984-FEB-5 1999-MAR-20 2008-APR-18 Altitude (ft) 4500 10400 21000 6800 Nominal Scale 9000 20000 40000 12000 Film Unknown Kodak 2405 DblX Kodak 2405 DblX Kodak 2405 DblX Camera  Unknown WILD RC-10 WILD RC-30 Wild RC-20 Focal Length (mm) Unknown 153.02 153.36 152.82 Lens Type Unknown 15 UAgII 15/4 UAG-S 15/4 UAG-S Lens Filter Unknown 525 nm (yellow) 525 nm (yellow) 525 nm (yellow) Number of Photos 35 18 13 48  All of the air photos scanned for this analysis include a margin, commonly called a titling bar, that contains useful information about the aerial survey the photos belongs to. Unfortunately, these margins confound subsequent image processing efforts if they are not removed. With the large number of photos included in this study, an automated image processing approach was developed to crop each photo by the interior dimensions of the photo margins (Appendix B.1).  After the images had been cropped, they were orthomosaiced using Agisoft Photoscan version 1.4.4 (Agisoft LLC, 2018), which is a digital photogrammetry software that relies on   17 structure from motion (SfM) to generate orthomosaics and 3D models of features or landscapes from a series of images that are captured with stereoscopic overlap. SfM software is particularly useful when analyzing historical imagery because it does not require information regarding original camera positioning (Fonstad, Dietrich, Courville, Jensen, & Carbonneau, 2013). The product of these efforts is a single mosaiced image object for each aerial survey (Appendix B.2).  It was a challenge to georeference the mosaics because reliable ground control points are unavailable for the four survey dates. Consequently, using ArcGIS, the 2008 mosaic was georeferenced to contemporary digital aerial images and the remaining mosaics were then sequentially georeferenced to one another in chronological order (Figure 4.7).    Figure 4.7 The workflow used to georeference each orthomosaic derived from archived aerial photographs.  4.3.2 Resample  The normal unit of measurement for spatial resolution in remote sensing is the ground sampled distance (GSD), which represents the smallest feature that can be reliably identified in an image (T. Nelson, Wulder, & Niemann, 2001; Thomson, 1994). True GSD of a scanned archived aerial photograph is determined and affected by the film, lens, scale, exposure, atmospheric haze, development process, scanning resolution, and any additional computational manipulation (Thomson, 1994). However, after an air photo has been scanned, spatial resolution is dictated only by scanning resolution and photo scale, and thus pixel size is not related to true GSD (T. Nelson et al., 2001).  Following the practical approach provided by Thomson (1994), the actual GSD of each aerial survey was measured by sampling vectors of intensity values perpendicular to 25 or more edges throughout each mosaic. Roof tops with dark shadows beneath them were preferentially chosen as sample locations given the stark contrast between high reflectance roof materials and low reflectance shadows. The elements within these vectors represent the intensity value of one pixel along each sample and together form an edge spread function (ESF), which captures the   18 amount of edge blurring at each sample location (Figure 4.8). Once these vectors were collected, the derivative between adjacent elements is computed to produce a line spread function (LSF) for each sample. The maximum slope value within the LSF represents the edge and the distance in either direction to the first instance of zero slope is the approximate linear extent of edge blurring. In theory, the LSF should be symmetrical (Thomson, 1994). However, it was observed that there was more variation in intensity values within low reflectance regions, which often produces a slightly asymmetrical or skewed LSF (Figure 4.8). Following Thomson (1994), the distance between the first instance of zero slope on either side of the edge is measured in meters and multiplied by 0.9 to produce an exclusive estimate of blurring. These estimates are then averaged to approximate the true GSD. Each mosaic was upscaled to match the calculated GSD using the bilinear interpolation function in ArcGIS (Table 4.2).    Figure 4.8 An example of a sample collected perpendicular to an edge (top) the resulting edge spread function (bottom-left) and the line spread function (bottom-right) produced by taking the derivative of the edge spread function.   19 Table 4.2 The ground sampling distance of each orthomosaic before and after they were resampled along with the nominal scale of each aerial survey  1951 1984 1999 2008 Nominal Scale 1:9000 1:20000 1:40000 1:12000 GSD Before Resample (m) 0.2463 0.4578 0.7432 0.2616 True GSD (m)  1.348  1.572 2.536 0.8703  4.3.3 Segmentation After resampling, the mosaics were masked by the outline of the active river channel so that ensuing classification efforts would not be confounded by features outside of the river. The example-based feature extraction function, provided by ENVI version 5.5 (Harris Geospatial Solutions, 2018), was then used to segment the resampled contemporary orthomosaic into sediment and non-sediment areas. This is an object-based approach, where objects or segments are formed by grouping pixels together based on their spectral, spatial, and/or texture attributes. This is in contrast with traditional pixel-based classification methods, where spectral information in each pixel is used to classify imagery. In theory, the approach can be fully automated with ground-truthed data, but in this case, due to the scale and simplicity of the classification and absence of such data, examples were provided manually based on visual assessment. After the mosaics were classified into their sediment and non-sediment components, the non-sediment regions were masked so that sediment regions could be analyzed on their own.   4.4 Contemporary Aerial Images 4.4.1 Acquisition, Orthomosaic, and Segmentation A set of contemporary digital aerial images with 18cm spatial resolution was captured by Dudley Thompson Mapping in early March 2018, during the Fraser River's low flow to ensure that the sediment bars present within the river were exposed to their greatest areal extent. These images were acquired using a Vexcel Ultracam Eagle large format digital camera mounted on a Piper Aztec PA23-250 twin engine turbo aircraft. The entire operation was confined to a single   20 12.9km flight line that yielded 12 images with 60% overlap and four spectral bands: red, green, blue, and near-infrared (Figure 4.9). Once the images were acquired they were orthomosaiced and georeferenced using Agisoft Photoscan to create a single object for further analysis (Appendix B.2). After the contemporary orthomosaic was constructed, areas inside the active channel were segmented into their sediment and non-sediment components using the example-based feature extraction function in ENVI.    Figure 4.9 The red line represents the ~13km flight used to collect 12 digital aerial images in March 2018. The blue circles denote the approximate principle point of each image while the yellow boundary represents the areal extent of the orthomosaic derived from the images.  4.4.2 Degradation The sediment regions of the contemporary orthomosaic were degraded so that they could become facsimiles of the sediment regions within the orthomosaics derived from each set of archived aerial photographs. While there are some similarities, notable differences in scale and exposure require degradation steps specific to each aerial survey. In an initial step, the multispectral contemporary orthomosaic was converted to panchromatic so that it would mimic the spectral characteristics of the archived aerial photographs. Aside from 1951, all of the aerial surveys used Kodak Double-X 2405 aerographic   21 film with a 525nm filter (Table 4.1). By filtering wavelengths up to 525nm, the film used during these surveys is sensitive only to wavelengths from 525 to ~715nm (Figure 4.10).   Figure 4.10 The spectral sensitivity of Kodak Double-X 2405 aerographic film (from Eastman Kodak Company, 2005), with an orange dashed line representing the 525 nm filter used in the 1984, 1999, and 2008 aerial surveys.  The sensor that measures red reflectance in the Vexcel Ultracam Eagle, which was used to capture the contemporary digital aerial images, is most sensitive between ~590 and 700 nm while the sensor that measures green reflectance is most sensitive between ~500 and 590 nm (Figure 4.11). The red band, therefore, spectrally covers ~65% of the sensitive range of the Kodak Double-X film and the green band covers ~35%. The following equation was used to convert the multispectral contemorary orthomosaic to panchromatic:  9:;&' = <=&' × 0.65A + <B&' × 0.35A  where i and j denote the column and row position of a pixel within the multispectral orthomosaic and the corresponding panchromatic orthomosaic. Unfortunately, all information concerning the film, camera, and lens used during the 1951 aerial survey was lost. As a consequence, the same equation used to mimic the spectral sensitivity of the other surveys was used to mimic the 1951 survey.   22  Figure 4.11 The spectral sensitivity of Vexcel Ultracam Eagle large format digital camera (from Vexcel Imaging, 2015), with an orange dashed line representing the 525 nm filter used in the 1984, 1999, and 2008 aerial surveys.   The bilinear interpolation function was then used within ArcGIS to alter the spatial resolution of the panchromatic contemporary orthomosaic so that it matched the GSD exhibited by each aerial survey. This produced four panchromatic orthomosaics each matching the spatial resolution of one of the historical aerial surveys. Finally, to mimic differences in exposure, the histogram of each degraded contemporary orthomosaic was matched to the histogram of the corresponding orthomosaic derived from archived aerial photographs in ArcGIS (Figure 4.12).     23  Figure 4.12 A schematic of the histogram matching transformation.  Following the recommendations of Geigel and Musgrave (1997) concerning the process of mimicking analog photographs, methods were explored to introduce granularity to the degraded contemporary images, which ultimately proved to be beyond the scope of this thesis. A lengthier discussion of grain simulation and an outline of possible future work will be provided in the discussion section.  4.5 Classification 4.5.1 Manual Interpretation A single contact print from each aerial survey was manually classified, by two skilled photo interpreters, into seven classes: gravel, sand, gravel with sand veneer, gravel/vegetation, sand/vegetation, water, and vegetation. The interpretations were recorded on semi-transparent plastic overlays and scanned using an EPSON Expression 120000XL scanner. The scans were digitized in ArcGIS so that they could be compared to one another as well as the computation-based classifications. Following the categorization outlined in Figure 4.13, two classifications were produced: one that delineated the sediment and non-sediment areas of each photo and one that outlined the gravel and sand within each photo.    24   Figure 4.13 The scheme used to categorize the manually interpreted classes.  4.5.2 Object-based The sediment regions of the degraded contemporary orthomosaics were classified into their gravel and sand components using an automated object-based image analysis (OBIA) routine. This routine involves grouping pixels into segments and then classifying those segments using georeferenced data based on their spectral, spatial, and textural attributes. Segments were delineated in ENVI using a watershed algorithm, which treats an image like a topographic map by defining “catchment basins” comprised of “ridges” of pixels with high intensity values and “troughs” of pixels with low intensity values. The size of these basins is determined by a scale parameter that ranges from 0 to 100 and controls the amount of “flooding” or maximum pixel intensity value that the “ridges” must contain. The result is a segmentation image, where each region contains the following attributes: mean intensity value, standard deviation of intensity values, size, compactness, and rectangularity. The degree of allowable difference between segments, based on distance and mean intensity value, can be controlled using a merge parameter which also ranges from 0 to 100.    25 The optimal combination of merge and scale parameters was unknown a priori and therefore the sediment regions of the degraded contemporary orthomosaics were segmented using each combination of merge and scale parameter at intervals of 10. This produced 121 segmented orthomosaics for each survey year. To determine the optimal parameter combination, the segmented orthomosaics were classified in ArcGIS using two robust, non-parametric, classifiers: support vector machine (SVM) and random tree (RTC). The classifications were trained and validated using stratified random samples comprised of 70% and 30% of the ground-truthed data, respectively. The optimal merge and scale parameters were chosen based on the combination that achieved the highest shared overall accuracy between the two classifiers. Once this occurred, the orthomosaic segmented with the optimal merge and scale parameters was classified using 100 bootstrapped training samples to determine the consistency of the two classifiers and the ground-truthed data set. In addition to classified imagery, these efforts produced a definition file for each bootstrapped sample, which parametrizes the corresponding classification. These definition files were used to classify scanned versions of the same archived aerial photographs that were manually interpreted.   4.5.3 Pixel-based In addition to object-based methods, a traditional pixel-based classification routine was also executed in ArcGIS based solely on spectral attributes. Again, SVM and RTC were used as classifiers. However, for this routine, a five meter buffer was computed around each ground-truthed location and stratified random samples of the pixels within those buffers were used as training and calibration data. The consistency of the classifiers was assessed using 100 bootstrapped samples. Similar to the object-based routines, the pixel-based classifiers produce a definition file which was used to classify the scanned contact prints that were manually interpreted.  4.5.4 Semi-automated  A semi-automated OBIA routine similar to the one used to segment the orthomosaics into their sediment and non-sediment components was employed to classify the sediment areas of each contact print into their sand and gravel components. The routine was carried out in ENVI   26 using the example-based feature extraction function, which segments the images using a watershed algorithm and then requires a user to provide examples of features to be classified. After the examples are provided, the segments are classified using a SVM classifier based on the attributes outlined in Table 4.3.   Table 4.3 The list of spectral, textural, and spatial attributes that are evaluated by the example-based feature extraction function in ENVI. Spectral Textural Spatial Mean Range Area Standard Deviation Mean Length Minimum Variance Compactness Maximum Entropy Convexity   Solidity   Roundness   Form Factor   Elongation   Rectangular Fit   Main Direction   Major Length   Minor Length   Number of Holes   Hole Area  The list of attributes available in this routine is much longer than the list of attributes that can be assessed in the automated object-based routine and more closely resembles the many characteristics that are evaluated during manual interpretation. Unfortunately, this routine does not produce a definition file that parameterizes the classification for reuse, and therefore requires that a new classification is developed for each new image or orthomosaic. Hence the moniker “semi-automated”.   4.5.5 Agreement The agreement between each of the manual interpretations as well as the automated and semi-automated classifications was assessed using a figure of merit defined by the following equation:   27  DEF	 = 	 94G94G	+	94H  where 94G is the number of mutually classified pixels and 94H is the number of uniquely  classified pixels (errors of commission and omission).     28 5 Results 5.1 Ground-truth Model The candidate models (Section 4.2.2) with the lowest Akaike information criterion (AIC) value at each unique number of degrees of freedom were selected and compared using an analysis of variance (ANOVA) test in R, which sequentially compares each model to determine if there is a significant reduction in residual deviance as independent variables are added. For example, the model with one independent variable is compared to the null model, the model with two independent variables is compared to the model with one independent variable, and so on. The results of this test, outlined in Table 5.1, show that the introduction of each new independent variable leads to a significant reduction in residual deviance at a = 0.01. Therefore, the most complex model, consisting of correlation, entropy, and brightness as fixed independent variables, was chosen for further analysis.   Table 5.1 The output from the analysis of variance test. Each model is compared to the next most complex model to determine if there is a significant reduction in residual deviance when an additional independent variable is added to the model.  Selected Models Residual Deviance Chi  Squared p value ~ (1|Day) -153.18904 NA NA ~ (1|Day) + Entropy -302.56588 149.376845 2.37E-34 ~ (1|Day) + Entropy + Correlation -320.97099 18.4051134 1.79E-05 ~ (1|Day) + Entropy + Correlation + Brightness -351.00181 30.0308192 4.25E-08  The predictions produced by the chosen model are plotted against observed percent gravel of the 120 training images in Figure 5.1. There appears to be good model fit throughout the range of observed values. By comparing a locally weighted smooth line with a 1:1 line, it becomes evident that there is a small amount of bias towards the middle and end of the range of observed values. The root mean square error (RMSE) of the model is 0.068 or 6.8%.    29  Figure 5.1 The predictions of percent gravel produced by the chosen model plotted against observed percent gravel. For comparison, a 1:1 line is provided (dotted orange) along with a locally weighted smooth line (solid blue).  The chosen model was bootstrap aggregated to improve accuracy, reduce variance, and avoid overfitting. The type of aggregator and the number of bootstrapped samples needed to outperform a single strong model was examined by computing the RMSE after each new sample had been aggregated. Median outperformed mean as an aggregator, and RMSE stabilized around 0.055 after ~750-1000 samples (Figure 5.2).     30  Figure 5.2 A comparison of mean and median as aggregators and the number of samples required to stabilize root mean square error.   5.2 Sediment and Non-sediment Agreement A semi-automated object-based classification routine was used to classify images into their sediment and non-sediment regions. The agreement between this routine and the two photo interpreters was assessed using a scanned contact print from each aerial survey. The two photo interpreters agree with one another at a slightly higher percentage than the semi-automated routine, but all agreement measures are above 90%. There does not seem to be any appreciable relationship between agreement measures and nominal scale, GSD, or survey year.            31 Table 5.2 The agreement between the semi-automated routine and the two photo interpreters when classifying an example archived aerial photograph from each survey into sediment and non-sediment.  1951  1984  Interp. #1 Interp. #2 Semi-automated  Interp. #1 Interp. #2 Semi-automated Interp. #1  0.97967686 0.95197215   0.96839392 0.90721715 Interp. #2   0.93654422    0.90183218         1999  2008  Interp. #1 Interp. #2 Semi-automated  Interp. #1 Interp. #2 Semi-automated Interp. #1  0.96729084 0.90950878   0.99959409 0.91862521 Interp. #2   0.9063318    0.91883723  5.3 Classification of Degraded Contemporary Orthomosaics The optimal scale and merge parameters, needed to segment the degraded contemporary orthomosaics (Section 4.5.2), are presented in Table 5.3. These are the parameters that produced the highest overall accuracy between two classifiers—support vector machine (SVM) and random tree (RTC). There does not appear to be a relationship between nominal scale or ground sample distance (GSD) and the scale and merge parameters. More comprehensive results depicting the overall accuracy and kappa coefficient for each merge and scale combination indicate that most combinations, where scale is between 0 – 50 and merge is between 0 – 80, will perform well (Appendix B.3). Beyond those ranges of merge and scale, performance drops off quickly as segments become too large.  Table 5.3 The optimal merge and scale parameters for each aerial survey. Year Scale Merge Ground Sampled Distance (m) Nominal Scale 1951 30 0 1.348 1:9000 1984 10 0 1.572 1:20000 1999 0 0 2.536 1:40000 2008 0 80 0.8703 1:12000   32 Once the optimal scale and merge parameters were identified, the two classifiers were bootstrapped using 100 training samples to assess their consistency. The overall accuracy of each of the bootstrapped classifications are presented in Figure 5.3. With the exception of a few outliers almost all classifications produced by the two classifiers have over 90% overall accuracy. These measures appear to be somewhat related to nominal scale with the 1951 degraded contemporary orthomosaic producing the best results and the 1999 degraded contemporary orthomosaic producing the worst results. Regardless of year, there does not seem to be any noticeable difference between the performance of the two classifiers.    Figure 5.3 The overall accuracy of 100 bootstrapped objected-based classifications for each degraded contemporary orthomosaic.   A similar process of bootstrapping SVM and RTC was used to assess their performance as pixel-based classifiers. The results of these classifications are provided in Figure 5.4. While the average overall accuracy of these classifications is slightly lower than the object-based classifications, there is a considerable level of consistency in performance across all survey   33 years. Again, similar to the object-based classifications there is very little difference in performance between the two classifiers.    Figure 5.4 The overall accuracy of 100 bootstrapped pixel-based classifications for each degraded contemporary orthomosaic.  5.4 Gravel and Sand Agreement The sediment regions within a single contact print from each aerial survey were manually classified into sand and gravel by two skilled photo interpreters. The same contact prints were analyzed using automated object- and pixel-based routines as well as a semi-automated object-based routine. The measures of agreement between these classifications, assessed using a figure of merit, are outlined in Table 5.4.  Agreement between the two photo interpreters appears to be related to survey year instead of scale with the highest measures of agreement produced for 2008 and the lowest for 1951. The automated pixel- and object-based routines have a high degree of agreement with one another (~85-98%) regardless of survey year, scale, or GSD, but low agreement with the manual interpretations (~45-70%). The semi-automated routine produced classifications with the most   34 balanced agreement between all of the different approaches. Within the aerial photographs from the 1999 and 2008 surveys, the semi-automated routine agreed with all of the different classifications between ~75-85%. For the 1951 photograph, the semi-automated routine produced low agreement with the manual interpretations (52% and 67%), but high agreement with the automated routines (~84-89%). In the 1984 photograph, the semi-automated routine produced measures of agreement between ~75-78% when compared to the manual interpretations and automated pixel-based classifications, but 71% and 67% when compared to the two automated object-based classifications.     35 Table 5.4 Agreement between the two photo interpreters, semi-automated routine, and automated pixel- and object-based routines when classifying an example archived aerial photograph from each survey into sand and gravel.  1951  Interp. #1 Interp. #2 Semi-automated Pixel - RTC Pixel - SVM OBIA - RTC OBIA - SVM Interp. #1  0.69350859 0.51875155 0.48139676 0.49768057 0.45345833 0.49294974 Interp. #2   0.6694524 0.66274469 0.67684139 0.62839421 0.65221928 Semi-automated    0.88224073 0.88982745 0.84170053 0.85977058 Pixel - RTC     0.97543834 0.8743421 0.8668212 Pixel - SVM      0.87739252 0.87516673 OBIA - RTC       0.9262772 OBIA - SVM        1:9000 Nominal Scale and 1.348m GSD  1984  Interp. #1 Interp. #2 Semi-automated Pixel - RTC Pixel - SVM OBIA - RTC OBIA - SVM Interp. #1  0.84292984 0.77528106 0.56417186 0.57726384 0.5391764 0.50949007 Interp. #2   0.77277398 0.54961229 0.56305058 0.52043175 0.47984521 Semi-automated    0.74824873 0.76300648 0.70614579 0.66980695 Pixel - RTC     0.98524221 0.87403937 0.86146614 Pixel - SVM      0.87144201 0.85856527 OBIA - RTC       0.93456612 OBIA - SVM        1:20000 Nominal Scale and 1.572m GSD     36  1999  Interp. #1 Interp. #2 Semi-automated Pixel - RTC Pixel - SVM OBIA - RTC OBIA - SVM Interp. #1  0.90937882 0.82141985 0.69187279 0.69891309 0.68608191 0.67293763 Interp. #2   0.86580755 0.72694119 0.73542203 0.71866421 0.70185201 Semi-automated    0.8254292 0.82981398 0.80464275 0.79778673 Pixel - RTC     0.97641495 0.85832828 0.85713421 Pixel - SVM      0.85685934 0.85532483 OBIA - RTC       0.96736676 OBIA - SVM        1:40000 Nominal Scale and 2.536m GSD   2008  Interp. #1 Interp. #2 Semi-automated Pixel - RTC Pixel - SVM OBIA - RTC OBIA - SVM Interp. #1  0.95662953 0.84398 0.68178072 0.71502908 0.63190865 0.62491775 Interp. #2   0.85186206 0.69120109 0.72495361 0.64456857 0.63743842 Semi-automated    0.81820558 0.85755736 0.76788571 0.75994726 Pixel - RTC     0.93430972 0.84256251 0.84777044 Pixel - SVM      0.85091284 0.85514187 OBIA - RTC       0.96930528 OBIA - SVM        1:12000 Nominal Scale and 0.8703m GSD     37 All of the bootstrapped pixel- and object-based definition files, derived from the degraded contemporary orthomosaics, were used to classify the contact prints to examine the consistency of the automated routines (Figure 5.5). Both routines produced stable average agreement measures at or above ~95% regardless of classifier, survey year, nominal scale, or GSD. Although differences are minimal, the pixel-based classifications appear to have slightly higher agreement along with greater consistency.  Figure 5.5 The agreement between 100 bootstrapped object-based (top) and pixel-based (bottom) classifications when applied to the sediment regions of an example archived aerial photograph from each survey   38 6 Discussion 6.1 Ground-truth Model The predictive model produced to determine the percent gravel within ground-level images containing a mixture of sand and gravel (Section 4.2.2) performed extremely well for the purposes of this study. The root mean square error (RMSE) for the model was 0.068 or 6.8%, which fell to ~0.055 or 5.5% after the model was incorporated into a bootstrap aggregating framework. Carbonneau et al. (2005) provide methods for accomplishing a similar objective in which the pixels of ground-level images are classified as either sand or gravel (Section 2.2). When replicated using the ground-level images collected for this study, their methods failed to produce a consistent result, largely due to sources of systematic variance such as the inconsistent angle and height at which the images were taken, variation in lighting conditions, and occasional presence of water or moisture on the bar surface. Object-based methods, outlined in Section 2.2, that are predicated on the work of Graham, Rice, & Reid (2005), were also explored by means of BASEGRAIN version 2.2 (Detert & Weitbrecht, 2012). Again, these methods failed to produce a consistent, reproducible results due to the multiple sources of systematic variation between images. This indicates that the model developed herein allows more flexibility during image capture than previously established methodologies. With that said, there is little doubt that removing one or more sources of systematic variance would improve model performance.   6.2 Degradation of Contemporary Digital Aerial Imagery The sediment regions of a contemporary multispectral orthomosaic, derived from digital aerial images, were degraded following the recommendations of Geigel & Musgrave (1997), so that they could be analyzed as a facsimile of the sediment regions within orthomosaics derived from archived aerial photographs (Section 4.4.2). The degradation process required that the contemporary orthomosaic was converted to panchromatic, resampled to match the ground sampling distance (GSD) of the historical aerial surveys, and transformed using a histogram matching routine. The underlying assumption of this routine is that the same variation in surficial sediment is present regardless of the year, and that the various characteristics of sediment such as reflectance and texture are relative. For example, if in one year gravel reflects highly in the visible spectrum relative to sand, it is assumed the same would be true any other year.    39 While the theoretical aspects of this approach follow established protocol, it is, unfortunately, nearly impossible to empirically test the efficacy of this process with the images available in this study. There are measures such as the structural similarity (SSIM) index that are designed to compare the similarity between an image and a degraded version of that image (Wang, Bovik, Sheikh, & Simoncelli, 2004). However, given that the archived aerial photographs and the degraded contemporary digital aerial images were captured decades apart, it would be difficult to determine how much of the difference between the two sets of images is due to flaws in the degradation process and how much is due to inherent differences in the phenomena captured in each set of images.  One additional shortcoming of the degradation process employed in this study is the exclusion of granularity. The presence of clumps of silver halide crystals within the emulsion of analog film results in a distinctive graininess, which should be simulated and introduced, if possible, to the digital images subject to degradation (Geigel & Musgrave, 1997). This is especially true for remote sensing efforts in which images are classified based on their spectral and textural attributes. Numerous quantitative methods have emerged to emulate the graininess of analog photographs in digital images. These methods range in complexity from white gaussian noise to robust, computationally expensive, Monte Carlo-based algorithms (A Newson, Delon, & Galerne, 2017). Unfortunately, these methods have been introduced without remote sensing products in mind, and assume that one pixel is comprised of multiple clumps of silver halide crystals. In the case of scanned analog aerial photographs, one clump of silver halide crystals can be represented by multiple pixels. The tremendous differences in the scale of hand-held images and remotely sensed products makes it difficult to easily employ established methods.  With multiple pixels representing individual clumps of silver halide crystals it was postulated that granularity could be quantified by measuring the degree of spatial autocorrelation between pixels using a 2D semivariogram. Once the semivariogram is generated, it would be incorporated into a geostatistical model to generate spatially correlated random fields (Dolloff & Doucette, 2014) that would then be overlaid, as a grain surface, multiplicatively or additively onto the contemporary digital images. While it is relatively straight-forward to generate semivarigorams for samples from each historical aerial survey, it became clear that justifying the   40 appropriate parameters values required for the geostatistical model and the methods to overlay the grain surface would be more difficult than justifying the absence of granularity.   6.3 Classifications of Degraded Contemporary Digital Aerial Imagery After the contemporary orthomosaic was degraded, automated pixel- and object-based classifications were trained and validated using 100 bootstrapped samples of ground-truthed data (Section 4.5.2 & 4.5.3). The object-based routines that evaluated spectral, textural, and spatial attributes produced overall accuracy measures ~1-2%  higher than the pixel-based routines, which only evaluated spectral intensity. Given the closeness of the results, this indicates that, at spatial resolutions of ~0.85 -2.5m, the defining characteristic of sand and gravel, as broad sedimentary classes, is spectral intensity instead of texture. An example of this phenomenon can be observed in Figure 6.1 where higher reflectance areas are indicative of gravel and lower reflectance areas are indicative of sand.   Figure 6.1 Points with ground-truthed estimates of percent gravel along with the contemporary orthomosaic after its degradation to match the aerial photographs from 2008. Each point represents a “site”, which is the mean center of four samples of percent gravel taken within a visually homogenous patch of sediment. The red label at each site location provides the average percent gravel value for the surrounding sediment patch.   41 6.4 Agreement Two experienced air photo interpreters manually classified one contact print from each historical aerial survey (Section 4.5.1). Both sets of four contact prints required approximately six hours of interpretation and an additional five hours to scan and digitize. Unsurprisingly, according to one of the interpreters, nominal scale played a large role in their ability to extract meaningful information from the prints regarding sand and gravel (Figure 6.2). True GSD calculated for each survey is an even better indication of the amount of detail that can be extracted from each contact print (Section 4.3.2). For example, while the 1951 survey was flown at the largest scale, the photos from that survey were printed on cardstock paper instead of photographic paper, resulting in a lower quality product and a level of detail similar to the 1984 survey.    Figure 6.2 The outline of the four contact prints used to assess agreement between various classification schemes.     42 All of the computational routines developed in this study to classify the archived aerial photographs were compared to the two manual interpretations to test their agreement (Section 4.5.5). The manual interpretations should be thought of as a standard instead of a reference. They represent what would be used to classify the archived aerial photographs if computational routines are not employed. It is therefore desirable if the computational routines similarly classify the same general areas, but define more accurate and distinct boundaries for each class.   6.4.1 Sediment and Non-sediment Classification The two interpreters produced preliminary classifications of sediment/non-sediment regions that achieved high levels of agreement (~97-100%) for all survey years. While slightly lower, the proposed semi-automated object-based routine still produced measures of agreement greater than 90% for all survey years. As expected, this indicates that, regardless of scale or GSD, sediment and non-sediment regions can be reliably classified based on visual assessment. The largest sources of discrepancy between manual interpretation and the proposed semi-automated routine appear to be the classification of transitional areas and boundary detail, both of which are exemplified in Figure 6.3. Visually, it is apparent that the semi-automated routine is far more detailed when defining boundaries than the manual interpretations. In transitional areas, small patches of sediment between larger areas of vegetation are classified as non-sediment by the manual interpreters and sediment by the semi-automated routine. Additionally, the semi-automated routine appears to mask out small pools of water that are difficult for the manual interpreters to distinguish from surrounding sediment areas.    43 Figure 6.3 Examples of sediment/non-sediment classification: original (top), manual interpretation (bottom left), semi-automated routine (bottom right). Non-sediment regions are in black.   44 6.4.2 Gravel and Sand Classification The two manual interpretations of sand and gravel produced measures of agreement that vary greatly depending on survey year. The contact print from 2008 was interpreted with ~95% agreement while the print from 1951 was interpreted with ~50% agreement. There does not seem to be a relationship between agreement and nominal scale or GSD. However, it is clear that at the level of detail provided by the manual interpretations, any divergence, even if it is a small, can have a significant impact on measures of agreement. Ultimately, lower agreement between manual interpretations decreases the reliability of any subsequent computational classification efforts for that survey year because it indicates that there is an absence of characteristics or sufficient detail that can be used to define areas of gravel and sand.   Figure 6.4 Examples of the two interpreter’s classification of the 2008 contact print into its sand and gravel components. Non-sediment regions are in black, sand is denoted by the colored regions, and gravel makes up the remaining areas of each image.   45 Figure 6.4 provides an example of the manual sand and gravel interpretations produced for the contact print from 2008. Both interpreters identify the same general areas as sand and gravel, but deviate in their class boundaries. Similar to the degraded contemporary orthomosaic, it can be observed that sand is represented by lower intensity values than gravel. This was captured by the automated computational routines, which also identified areas with lower reflectance as sand and areas with higher reflectance as gravel (Figure 6.5). Across all survey years, the automated routines agreed with one another between ~85-98%. There does not appear to be a noticeable difference between SVM or RTC as classifiers for either classification paradigm. When compared to the manual interpretations, agreement measures ranged from ~45% for the 1951 print to ~70% for the 1999 and 2008 prints. Agreement between manual and automated classification efforts seems to be, qualitatively, related to the level of agreement between manual interpreters. None of these results indicate that the automated object-based routines provide any appreciable benefits over the pixel-based routine, especially considering that additional computation is required to execute object-based classifications. Neither classification paradigm is noticeably more consistent, in greater agreement with the manual interpretations, or produced higher overall accuracies in the contemporary orthomosaic classifications (Section 5.3 & 5.4). This suggests that pixel-based routines should be employed for future efforts, if the methodology outlined in this study is replicated and automated routines are desired.    46  Figure 6.5 An example of each automated routine used in this study: object-based SVM classification (top left), object-based RTC classification (top right), pixel-based SVM classification (bottom left), pixel-based RTC classification (bottom right).   47  The final classification that was examined in this study is a semi-automated object-based routine that was carried out completely within ENVI version 5.5 using the example-based feature extraction function (Section 4.5.4). The agreement measures produced by this classification indicate that it is somewhat of a hybrid between manual interpretation and the automated classifications. This is likely because the classification process is supervised by a user  who can provide training examples that balance measurable spectral, textural, and spatial attributes along with more complex attributes such as site, pattern, shadow, and association.     The primary shortcoming of this routine is its reproducibility. The semi-automated classification must be retrained for each new mosaic or image while the automated classifications generate a definition file that can be reused to classify additional images or orthomosaics. As a result, the automated routines are not subject to the same inconsistencies due to human error. They also do not require an individual who is experienced with aerial photography and image processing. With that said, if they are available, an inexperienced user could reference a manually interpreted contact print or ground-truthed data indicating the locations of sand and gravel within a degraded contemporary orthomosaic in order to produce a Figure 6.6 An example of the semi-automated classification of the 2008 contact print into its sand and gravel components. Non-sediment regions are in black, sand is denoted in orange, and gravel makes up the remaining areas of the image.   48 classification that is potentially very similar in accuracy to a classification produced by a more experienced user.  6.5 Future Work 6.5.1 Ground-truth Model The parametric model that was developed in this study to generate estimates of percent gravel in ground-level images relies on the beta probability distribution (Section 4.2.2). Such models attempt to define a data distribution using a finite set of parameters. This makes them less flexible than non-parametric models, which allow for an infinite number of parameters that are dynamic enough to grow and change along with the training data. Non-parametric models such as neural networks, support vector machines, and random trees (forest) are particularly useful when predictions rather than inference are of interest. Now that the relationship between image texture and grain size has been established, in future studies, it may be prudent to abandon the long standing use of parametric models in favor of non-parametric algorithms.   6.5.2 Additional Variables for Classification In this study, automated and semiautomated pixel- and object-based routines were explored. The pixel-based routines relied solely on spectral information while the object-based routines relied on textural and spatial attributes in addition to spectral information. The list of attributes humans can consider when classifying an image is much longer than either of these computational classification paradigms (Figure 6.7). Additional attributes could be introduced to object- or pixel-based classifications that mimic the attributes evaluated during manual classifications. For example, in this study, attributes providing information regarding the orientation or position of a pixel or segment in relation to channel flow or the channel bar could be useful in their classification as well as a more theoretical understanding of sand and gravel deposition. Height derived from the 3D model used to orthorectify the contemporary and archived mosaics could also be used to more effectively classify images into their sediment and non-sediment components since there is an obvious height difference between vegetation, water, and sediment. The use of height could also be extended to the gravel and sand classification with some success, but the efficacy of this application is less assured.    49   Figure 6.7 The primary ordering of image analysis elements that were used by each classification paradigm in this study (adapted from Estes, Hajic, & Tinney, 1983).  6.5.3 Degradation 6.5.3.1 Granularity As stated earlier, granularity was not included in the degradation process developed in this study. More work is required to establish a consistent grain simulation methodology that can be applied to remotely sensed products at varying scales. This might involve the adaptation of established routines for images captured with hand-held cameras or it may require something new. A geostatistical approach in which a grain surface is simulated using a 2D semivariogram and overlaid onto digital images has some promise as a computationally inexpensive and scalable routine. However, the appropriate geostatistical model parameters along with the overlay process need additional examination.  6.5.3.2 Generative Adversarial Neural Networks Generative adversarial neural networks (GANs) are an exciting new class of machine learning algorithm that can be trained to learn the attributes of images, paintings, and music in order to reproduce those attributes to create facsimiles (Goodfellow et al., 2014). Recently, a GAN was used to learn the characteristics of paintings by various well-known painters and then   50 reproduce those learned characteristics in digital images (Tan, Chan, Aguirre, & Tanaka, 2017). In a sense, the degradation process carried out in this study was done “by hand” with each step parameterized and executed using code. Theoretically, it is feasible that this process could be automated if the attributes of archived aerial photographs are learned and then replicated by a GAN in digital aerial imagery.  6.5.3.3 Assessment For true, quantitative measures of efficacy, degradation processes would need to be assessed by capturing digital and analog images at roughly the same time. One could then use measures such as structural similarity after each successive degradation step to determine the importance of each step in the overall process. It is recommended that such work be carried out within a study area containing features and structures that are time invariant and relatively homogenous in their spectral and textural attributes.                              51 7 Conclusion This study provides an exploration of methods capable of delineating patches of sand and gravel within archived aerial photographs. Automated classification routines were developed by pairing ground-truthed data with contemporary digital aerial imagery that was degraded to mimic the attributes of the photos from each aerial survey. The automated classifications were compared to semi-automated object-based classifications as well as to manual interpretations of photos from each survey. The main conclusions are outlined below. • Ground-truth Model: Measures of image texture derived from a grey-level co-occurrence matrix are strongly related to the amount of surficial gravel within ground-level images. This relationship can be exploited to develop a predictive generalized linear model that performs well regardless of the consistency with which the ground-level images are captured.  • Degraded Contemporary Digital Aerial Imagery: o Patches of sand and gravel within contemporary digital aerial imagery that has been degraded to mimic the characteristics of analog aerial photography can be consistently classified with overall accuracy above ~93% using automated object- and pixel-based classifiers.  o At spatial resolutions of ~0.85-2.5m, the defining characteristic of sand and gravel patches is spectral intensity instead of image texture. • Scanned Archived Aerial Photographs: o Object-based classifications do not provide an appreciable increase in overall accuracy or consistency when used to delineate sand and gravel patches within contemporary digital aerial imagery or scanned analog photographs. Automated pixel-based classifications that rely solely on spectral information should be used instead.    52 o Sediment and non-sediment regions can be classified within archived aerial photographs using a semi-automated object-based routine. These classifications have ~90% agreement with manual interpretations.  o Semi-automated object-based classifications of sand and gravel within archived aerial photographs provide the best balance between manual interpretation and automated computational routines, but they must be retrained each time a new image or orthomosaic needs to be analyzed. This can result in inconsistency, especially if the classifications are trained by someone who is inexperienced with fluvial geomorphology and image processing.                  53 References Agisoft LLC. (2018). Agisoft PhotoScan Professional (Version 1.4.4). Retrieved from http://www.agisoft.com/downloads/installer/ Black, M., Carbonneau, P., Church, M., & Warburton, J. (2014). Mapping sub-pixel fluvial grain sizes with hyperspatial imagery. Sedimentology, 61(3), 691–711. https://doi.org/10.1111/sed.12072 Brasington, J., Langham, J., & Rumsby, B. (2003). Methodological sensitivity of morphometric estimates of coarse fluvial sediment transport. Geomorphology, 53(3–4), 299–316. https://doi.org/10.1016/S0169-555X(02)00320-3 Breiman, L. (1996). Bagging predictors. Machine Learning, 24(2), 123–140. https://doi.org/10.1007/BF00058655 Buffington, J. M., & Montgomery, D. R. (1999). A procedure for classifying textural facies in gravel-bed rivers. Water Resources Research, 35(6), 1903–1914. https://doi.org/10.1029/1999WR900041 Carbonneau, P. E., Bergeron, N. E., & Lane, S. N. (2005). Texture-based image segmentation applied to the quantification of superficial sand in salmonid river gravels. Earth Surface Processes and Landforms, 30(1), 121–127. https://doi.org/10.1002/esp.1140 Carbonneau, P. E., Lane, S. N., & Bergeron, N. E. (2004). Catchment-scale mapping of surface grain size in gravel bed rivers using airborne digital imagery. Water Resources Research, 40(7), 1–11. https://doi.org/10.1029/2003WR002759 Chandler, J. H., Rice, S., & Church, M. (2004). Colour aerial photography for riverbed classification. In 20th ISPRS Congress on Technical Commission VII (Vol. 35, pp. 1079–1084). Retrieved from https://www.scopus.com/inward/record.uri?eid=2-s2.0-  54 77954379752&partnerID=40&md5=001e04bf5e73de8c293f70885292b244 Church, M. (2006). Bed Material Transport and the Morphology of Alluvial River Channels. https://doi.org/10.1146/annurev.earth.33.092203.122721 Church, M. (2017). Landscapes and Landforms of Western Canada. In O. Slaymaker (Ed.), World Geomorphological Landscapes (pp. 381–393). Cham, Switzerland: Springer. https://doi.org/10.1007/978-3-319-44595-3 Clausi, D. A. (2002). An analysis of co-occurrence texture statistics as a function of grey level quantization. Canadian Journal of Remote Sensing, 28(1), 45–62. https://doi.org/10.5589/m02-004 Cliff, A. D. (Andrew D., & Ord, J. K. (1973). Spatial autocorrelation. Pion. Retrieved from https://catalogue.nla.gov.au/Record/608160 Cots-Folch, R., Aitkenhead, M., & Martinez-Casasnovas, J. (2007). Mapping land cover from detailed aerial photography data using textural and neural network analysis. International Journal of Remote Sensing, 28(7–8), 1625–1642. https://doi.org/10.1080/01431160600887722 Cullen, A. C., & Frey, C. H. (1999). Probabilistic techniques in exposure assessment: a handbook for dealing with variability and uncertainty in models and inputs. Springer Science & Business Media. Detert, M., & Weitbrecht, V. (2012). Automatic object detection to analyze the geometry of gravel grains - A free stand-alone tool. River Flow 2012 - Proceedings of the International Conference on Fluvial Hydraulics, 1, 595–600. Dolloff, J., & Doucette, P. (2014). The Sequential Generation of Gaussian Random Fields for Applications in the Geospatial Sciences. ISPRS International Journal of Geo-Information,   55 3(2), 817–852. https://doi.org/10.3390/ijgi3020817 Environment Canada. (2016). Monthly Discharge Graph for FRASER RIVER AT HOPE (08MF005) [BC]. Retrieved June 30, 2019, from https://wateroffice.ec.gc.ca/report/historical_e.html?stn=08MF005&mode=Graph&type=h2oArc&results_type=historical&dataType=Monthly&parameterType=Flow&year=2016&y1Max=1&y1Min=1 Environmental Systems Research Institute. (2016). ArcGIS Release 10.5. Redlands, CA. Fonstad, M. A., Dietrich, J. T., Courville, B. C., Jensen, J. L., & Carbonneau, P. E. (2013). Topographic structure from motion: A new development in photogrammetric measurement. Earth Surface Processes and Landforms, 38(4), 421–430. https://doi.org/10.1002/esp.3366 Garmin. (2011). eTrex ® Owner’s Manual. Retrieved from www.garmin.com Geigel, J., & Musgrave, F. K. (1997). A model for simulating the photographic development process on digital images. In Proceedings of the 24th annual conference on Computer graphics and interactive techniques  - SIGGRAPH ’97 (pp. 135–142). New York, New York, USA: ACM Press. https://doi.org/10.1145/258734.258813 Gonzalez, R. C., & Woods, R. E. (2008). Digital Image Processing. Prentice Hall (3rd ed.). Upper Saddle River, New Jersey: Pearson Prentice Hall. https://doi.org/10.1016/0734-189X(90)90171-Q Goodfellow, I. J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., … Bengio, Y. (2014). Generative Adversarial Nets. In Advances in Neural Information Processing Systems. Retrieved from http://www.github.com/goodfeli/adversarial Graham, D. J., Reid, I., & Rice, S. P. (2005). Automated sizing of coarse-grained sediments: Image-processing procedures. Mathematical Geology, 37(1), 1–28.   56 https://doi.org/10.1007/s11004-005-8745-x 10.1001/jama.288.23.3035 Graham, D. J., Rice, S. P., & Reid, I. (2005). A transferable method for the automated grain sizing of river gravels. Water Resources Research, 41(7), 1–12. https://doi.org/10.1029/2004WR003868 Ham, D. (2005). Morphodynamics and sediment transport in a wandering gravel-bed channel : Fraser River , British Columbia. University of British Columbia. Haralick, R. (1979). Statistical and structural approach to texture. Proceeding of IEEE Vol 67 No 5, 67(5), 786–804. https://doi.org/10.1109/PROC.1979.11328 Haralick, R., Shanmugan, K., & Dinstein, I. (1973). Textural features for image classification. IEEE Transactions on Systems, Man and Cybernetics. https://doi.org/10.1109/TSMC.1973.4309314 Harris Geospatial Solutions. (2018). ENVI (version 5.5). Ibbeken, H., & Schleyer, R. (1986). Photo-sieving: A method for grain-size analysis of coarse-grained, unconsolidated bedding surfaces. Earth Surface Processes and Landforms, 11(1), 59–77. https://doi.org/10.1002/esp.3290110108 Kelly, D. H. (1960). Systems Analysis of the Photographic Process I A Three-Stage Model. Journal of the Optical Society of America, 50(3), 269. https://doi.org/10.1364/josa.50.000269 Eastman Kodak Compay. (2005). Kodak Double-X Aerographic Film 2405. Lane, S. N., Widdison, P. E., Thomas, R. E., Ashworth, P. J., Best, J. L., Lunt, I. A., … Simpson, C. J. (2010). Quantification of braided river channel change using archival digital image analysis. Earth Surface Processes and Landforms, 35(8), 971–985. https://doi.org/10.1002/esp.2015   57 Marcus, W. A., & Fonstad, M. A. (2008). Optical remote mapping of rivers at sub-meter resolutions and watershed extents. Earth Surface Processes and Landforms, 33(1), 4–24. https://doi.org/10.1002/esp.1637 Morgan, J. L., & Gergel, S. E. (2010). Quantifying historic landscape heterogeneity from aerial photographs using object-based analysis. Landscape Ecology, 25(7), 985–998. https://doi.org/10.1007/s10980-010-9474-1 Morgan, J. L., Gergel, S. E., & Coops, N. C. (2010). Aerial Photography: A Rapidly Evolving Tool for Ecological Management. BioScience, 60(1), 47–59. https://doi.org/10.1525/bio.2010.60.1.9 Neill, C. R. (1987). Sediment Balance Considerations Linking Long-term Transport and Channel Processes. In C. R. Thorne, J. C. Bathurst, & R. D. Hey (Eds.), Sediment Transport in Gravel-Bed Rivers. Chichester: J. Wiley. Retrieved from https://trove.nla.gov.au/work/18946022?q&versionId=22239682 Nelson, P., Bellugi, D., & Dietrich, W. (2014). Delineation of river bed-surface patches by clustering high-resolution spatial grain size data. Geomorphology, 205, 102–119. https://doi.org/10.1016/j.geomorph.2012.06.008 Nelson, T., Wulder, M., & Niemann, O. (2001). Spatial resolution implications of digitizing aerial photography for environmental applications. The Imaging Science Journal, 49(4), 223–232. https://doi.org/10.1080/13682199.2001.11784386 Newson, A, Delon, J., & Galerne, B. (2017). A Stochastic Film Grain Model for Resolution-Independent Rendering. COMPUTER GRAPHICS Forum, 36(8), 684–699. https://doi.org/10.1111/cgf.13159 Newson, Alasdair, Faraj, N., Delon, J., & Galerne, B. (2017). Analysis of a Physically Realistic   58 Film Grain Model, and a Gaussian Film Grain Synthesis Algorithm. In 6th Conference on Scale Space and Variational Methods in Computer Vision. Kolding, Denmark. Retrieved from https://hal.archives-ouvertes.fr/hal-01494123v2 Okeke, F., & Karnieli, A. (2006a). Methods for fuzzy classification and accuracy assessment of historical aerial photographs for vegetation change analyses. Part I: Algorithm development. International Journal of Remote Sensing, 27(1), 153–176. https://doi.org/10.1080/01431160500166540 Okeke, F., & Karnieli, A. (2006b). Methods for fuzzy classification and accuracy assessment of historical aerial photographs for vegetation change analyses. Part II: Practical application. International Journal of Remote Sensing, 27(9), 1825–1838. https://doi.org/10.1080/01431160500315691 Rice, S., & Church, M. (1996). Sampling Surficial Fluvial Gravels: The Precision of Size Distrubution Precentile Estimates. Journal of Sedimentary Research, 66(3), 654–665. Rice, S., & Church, M. (2010). Grain-size sorting within river bars in relation to downstream fining along a wandering channel. Sedimentology, 57(1), 232–251. https://doi.org/10.1111/j.1365-3091.2009.01108.x Rubin, D. M. (2004). A Simple Autocorrelation Algorithm for Determining Grain Size from Digital Images of Sediment. Journal of Sedimentary Research, 74(1), 160–165. https://doi.org/10.1306/052203740160 Tan, W. R., Chan, C. S., Aguirre, H. E., & Tanaka, K. (2017). ArtGAN: Artwork synthesis with conditional categorical GANs. In 2017 IEEE International Conference on Image Processing (ICIP) (pp. 3760–3764). IEEE. https://doi.org/10.1109/ICIP.2017.8296985 Thomson, G. H. (1994). A Practical Method of Determining the Ground Sampled Distance in   59 Small Scale Aerospace Photography. The Journal of Photographic Science, 42(2), 64–64. https://doi.org/10.1080/00223638.1994.11738564 Venditti, J. G., & Church, M. (2014). Morphology and controls on the position of a gravel- sand transition: Fraser River, British Columbia. Journal of Geophysical Research: Earth Surface, 119(9), 1959–1976. Verdú, J. M., Batalla, R. J., & Martínez-Casasnovas, J. A. (2005). High-resolution grain-size characterisation of gravel bars using imagery analysis and geo-statistics. Geomorphology, 72(1–4), 73–93. https://doi.org/10.1016/j.geomorph.2005.04.015 Vexcel Imaging. (2015). UltraCam Eagle Calibration Report. Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2004). Image Quality Assessment: From Error Visibility to Structural Similarity. IEEE TRANSACTIONS ON IMAGE PROCESSING (Vol. 13). Retrieved from http://www.cns.nyu.edu/~lcv/ssim/. Wentworth, C. K. . (1922). A Scale of Grade and Class Terms for Clastic Sediments. The Journal of Geology, 30(5), 377–392. Wolman, M. G. (1954). A method of sampling coarse river-bed material. Eos, Transactions American Geophysical Union, 35(6), 951–956. https://doi.org/10.1029/TR035i006p00951          60 Appendices Appendix A    This Appendix outlines complementary material for the methods used to analyze and process ground-level images collected in October 2017 and March 2018.  A.1 Ground-level Image Preprocessing All of the ground-level images of surficial sediment contained a 0.75m X 1.0m wooden frame. An automated image processing routine was developed in MATLAB so that each image could be cropped by the inner dimensions of the wooden frame. This ensured that all of the images represented the same approximate areal extent and that the frame wouldn’t confound any subsequent analysis. A visual representation of this routine is provided in Figure A.1 and the code is given on the following page.      Figure A.1 A visual representation of the image processing routine used to crop each of the ground-level images.    61 imageDir = Input folder directory  TestDir = Output folder directory filePattern = fullfile(imageDir, '*.JPG'); ImageFiles = dir(filePattern);  for k = 1:length(ImageFiles) %Read image and establish output protocol baseFileName = ImageFiles(k).name; fullFileName = fullfile(imageDir, baseFileName); TestSave = fullfile(TestDir, baseFileName); I = imread(fullFileName);      %Isolate the frame and then crop by its interior dimensions  Gry = rgb2gray(I); %rgb to greyscale image contrast = adapthisteq(Gry); %histogram stretch to increase constrast se = strel('square', 10); erodedI = imerode(contrast, se); %erode image level = graythresh(erodedI);  BW = imbinarize(erodedI, level); %binarize eroded image bw = bwareafilt(BW, 1); %select the large object (should be frame) invert = imcomplement(bw); %invert J_area = bwareafilt(invert, 1); %select the large object (should be are within the frame)      %Get bounding image coordinates of a convex hull that encompasses the area %within the frame. Crop by those coordinates bw_convex = regionprops(J_area, 'ConvexHull');  convex = cell2mat(struct2cell(bw_convex)); [~, ix,] = min(sum(convex, 2)); xmin_ymin = convex(ix, 1:2); xmin = xmin_ymin(:,1); ymin = xmin_ymin(:,2); [CM, ~] = max(convex(:,1)); [RM, ~] = max(convex(:,2)); width = CM - xmin; height = RM - ymin; rect = [xmin, ymin, width, height]; croppedI = imcrop(I, rect);      %Now that the image has been cropped by the frame we know the y-axis is ~1 meter in length, which allows us to compute the number of pixels per SI unit [x, y] = size(croppedI(:,:,1)); %Get dimensions of image PxCm = round(y/100); %Pixels per centimeter nail = PxCm * 6; %The nails are 6 cm long, so multiply PxCm by 6      %Get xmin, ymin and xmax,ymax based on the length of the nails. And crop ColMax = y - nail;  RowMax = x - nail; xmin = nail; ymin = nail; width = ColMax - xmin; height = RowMax - ymin;  rect = [xmin, ymin, width, height]; croppedI = imcrop(croppedI, rect); imwrite(croppedI, TestSave); end   62 A.2 Texture Statistics Generated from Grey-level Co-occurrence Matrix A bootstrap aggregated generalized linear model was developed to predict the proportion of gravel within ground-level images containing a mixture of sand and gravel. The independent variables for this model were computed for each image using MATLAB.  Five of the six independent variables calculated for each image are texture statistics derived from a grey-level co-occurrence matrix (GLCM). These texture statistics and their formulas are outlined in Table A.1  Table A.1 Texture statistics generated using a GLCM (adapted from Gonzalez & Woods, 2008) Statistics Description Formula Contrast Measures the local variations in the GLCM !!(i − j)'() p)( Correlation Measures the joint probability occurrence of the specified pixel pairs 	∑ ∑ (- − .)/(-.)01 − µ3µ4σ3σ4  Energy (Angular Second Moment or Uniformity) Provides the sum of squared elements in the GLCM !!(- − .)'01  Entropy A measure of randomness within the GLCM −!!/106789/10:01  Homogeneity Measures the closeness of the distribution of elements in the GLCM to the GLCM diagonal. !! /101 + |- − .|01              63 The MATLAB code used to compute the independent variables for each ground-level image:  FolderDir = %Input folder directory containing cropped ground-level images FolderOI = dir(FolderDir); FolderOI(17).name); FolderPattern = fullfile(FolderDir, '*.jpg'); ImageFiles = dir(FolderPattern); a = struct('Contrast', {}, 'Correlation', {}, 'Energy', {},'Homogeneity', {}, 'Entropy', {},'Name', {}, 'Red', {}, 'Green', {}, 'Blue', {}, 'Std', {}, 'Day', {});  for k = 3:length(ImageFiles)     baseFileName = ImageFiles(k).name;     fullFileName = fullfile(ImageFiles(k).folder, baseFileName);     I = imread(fullFileName);     new = rgb2gray(I);     glcm = graycomatrix(new);     stats = graycoprops(glcm);     stats.Name = baseFileName(1:end-4);     position = k - 2;     E = entropy(new);     [rows, cols, ~] = size(I);     band_means = sum(sum(I,2),1)/(rows*cols);      a(position).Name = stats.Name;     a(position).Contrast = stats.Contrast;     a(position).Correlation = stats.Correlation;     a(position).Energy = stats.Energy;     a(position).Homogeneity = stats.Homogeneity;     a(position).Entropy = E;     a(position).Red = band_means(1);     a(position).Blue = band_means(2);     a(position).Green = band_means(3);     a(position).Std = std2(new);     a(position).Day = baseFileName(1:2);      %Output GLCM stats as .csv     TablePath = fullfile('~\\GLCM_Stats\\', FolderOI(17).name);     TableExt = '.csv';     TableName = strcat(TablePath, TableExt);     writetable(struct2table(a), TableName); end                64 A.3 Generalized Linear Model The output for the best performing generalized linear model including parameters and metrics:  Family: beta  (logit) Formula:  PercentGravel ~ Entropy + Correlation + Brightness + (1 | Day) Data: train  AIC      BIC   logLik deviance df.resid -339.0   -322.1    175.5   -351.0      117  Conditional model: Groups Name        Variance Std.Dev. Day    (Intercept) 0.3386   0.5819 Number of obs: 120, groups:  Day, 6  Overdispersion parameter for beta family (): 38.4  Conditional model: Estimate Std. Error z value Pr(>|z|) (Intercept) -33.111816   1.937184 -17.093  < 2e-16 *** Entropy       4.337714   0.282633  15.348  < 2e-16 *** Correlation   4.263472   0.731150   5.831 5.50e-09 *** Brightness   -0.010913   0.001868  -5.842 5.15e-09 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’                      65 The R code that was written to develop and deploy the bootstrap aggregated generalized linear model:  library(ggplot2) library(lme4) library(DHARMa) library(fitdistrplus) library(boot) library(ResourceSelection) library(rcompanion) library(lmtest) library(car) library(MuMIn) library(dplyr) library(foreach) library(matrixStats)  setwd("F:\\Results\\GroundTruthedModel")  ####Read data  df <- read.csv("TextureStatistics.csv") df$Brightness <- as.numeric((df$Green + df$Red + df$Blue)/3) #Create brightness variable train <- df %>%   filter(PercentGravel > 0.01 & PercentGravel < 0.99) #Remove samples that contain 100% gravel/sand  ####Cullen and Frey Graph  descdist(train$PercentGravel, boot = 10000, method = 'sample')  ####Multiple Beta Regression Test library(glmmTMB) selecteddf <- train[,c(1:5,13)] #Select texture statistics and brightness as variables to test myVariables <- names(selecteddf) null <- "PercentGravel ~ (1|Day)" #Specify null model out <- unlist(lapply(1:5, function(n) combn(myVariables, n,                            FUN=function(row) paste0("PercentGravel ~ (1|Day)+",                                                     paste0(row, collapse = "+"))))) #Create vector of variable combinations up to combinations of five null[2:(length(out)+1)] <- out myFormulas <- as.list(null)  #Fit a model for each set of variable combinations myBetas <- list() for(i in 1:length(myFormulas)) {   myBetas[[i]] <- glmmTMB(as.formula(myFormulas[[i]]), family = beta_family(link = "logit"), data = train) } names(myBetas) <- null[1:length(myBetas)]  #Compute some statistics for each model   66 RMSE <- vector() Hoslem <- vector() aic <- vector() loglike <- vector() dev <- vector() df_resid <- vector() PseudoR2 <- vector() for(i in 1:length(myBetas)) {   train$BetaYhat <- fitted(myBetas[[i]])   train$BetaError <- train$PercentGravel - train$BetaYhat      RMSE[i] <- sqrt(mean((train$BetaError)^2))   Hoslem[i] <- hoslem.test(train$PercentGravel, train$BetaYhat)[["p.value"]]   aic[i] <- AIC(myBetas[[i]])   dev[i] <- deviance(myBetas[[i]])   df_resid[i] <- df.residual(myBetas[[i]])   loglike[i] <- logLik(myBetas[[i]])   PseudoR2[i] <- cor(train$PercentGravel,train$BetaYhat)^2 }  #Compine statistics into a table with model formula ModelComparison <- data.frame(cbind(null[1:length(myBetas)], PseudoR2, RMSE, Hoslem, aic, loglike, dev, df_resid)) ModelComparison$aic <- as.numeric(as.character(ModelComparison$aic)) ModelComparison$RowNames <- row.names(ModelComparison)  SmallestAIC <- ModelComparison %>%    group_by(df_resid) %>% summarise(Value = min(aic)) names(SmallestAIC)[2] <- "aic"  SmallestAIC <- left_join(SmallestAIC, ModelComparison, by = c("df_resid", "aic"))  comparison <- anova(myBetas[[as.numeric(SmallestAIC[1,"RowNames"])]], myBetas[[as.numeric(SmallestAIC[2,"RowNames"])]],        myBetas[[as.numeric(SmallestAIC[3,"RowNames"])]], myBetas[[as.numeric(SmallestAIC[4,"RowNames"])]],        myBetas[[as.numeric(SmallestAIC[5,"RowNames"])]], myBetas[[as.numeric(SmallestAIC[6,"RowNames"])]])  comparison$Model <- rev(SmallestAIC$V1) row.names(comparison) <- 1:nrow(comparison) comparison <- comparison[,c(9,1:8)] write.csv(comparison, "F:\\Results\\GroundTruthedModel\\ExploratoryComparison.csv")  anova(myBetas[[as.numeric(SmallestAIC[1,"RowNames"])]], myBetas[[as.numeric(SmallestAIC[3,"RowNames"])]])  ####Test For Colinearity library(mctest)  X <- train[,c(1:5,13)]  omcdiag(X, train$PercentGravel) imcdiag(X, train$PercentGravel)   67  ####Chosen Mixed Effects Beta Model library(glmmTMB)  GlmmBeta <- glmmTMB(PercentGravel ~ Entropy + Correlation + Brightness + (1|Day), data = train, family = beta_family(link = "logit")) summary(GlmmBeta)  save(GlmmBeta, file = "F:\\Results\\GroundTruthedModel\\GlmmBeta") capture.output(summary(GlmmBeta), file = "F:\\Results\\GroundTruthedModel\\GlmmBeta.txt")  train$BetaYhat <- fitted(GlmmBeta) train$BetaError <- train$PercentGravel - train$BetaYhat  sqrt(mean((train$BetaError)^2)) hoslem.test(train$PercentGravel, train$BetaYhat)  ggplot(data = train) +   geom_point(aes(PercentGravel, BetaYhat)) +   geom_smooth(aes(PercentGravel, BetaYhat), method = "lm") +    scale_y_continuous(breaks=c(0,0.25,0.5,0.75,1)) +   scale_x_continuous(breaks=c(0,0.25,0.5,0.75,1)) +   labs(title = "Predicted Percent Gravel in Ground Truthed Images",         y = "Predicted Percent Gravel", x = "Observed Percent Gravel")  ggplot(train, aes(BetaYhat, BetaError)) +   geom_point() +   geom_smooth(method = "lm") +   scale_x_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1)) +   scale_y_continuous(breaks=c(-.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +   labs(title = "Error Variance of Percent Gravel Predictions",         y = "Error", x = "Predicted Percent Gravel")  ####RMSE test df <- read.csv("TextureStatistics.csv") df$Brightness <- as.numeric((df$Green + df$Red + df$Blue)/3) train <- df %>%   filter(PercentGravel > 0.01 & PercentGravel < 0.99)  positions <- sample(nrow(train),size=floor((nrow(train)/4)*3)) training<- train[positions,] testing<- train[-positions,]  AvgRmse <- list() MedianRmse <- list() iterations<-2000 predictions <- foreach(m=1:iterations,.combine=cbind) %do%  {   training_positions <- sample(nrow(training), size=floor((nrow(training)/4)))   train_pos<-1:nrow(training) %in% training_positions   trainingdf <- training[train_pos,]      GlmmBeta <- glmmTMB(PercentGravel ~ Entropy + Brightness + Correlation + (1|Day), data = trainingdf, family = beta_family(link = "logit"))      68   predict(GlmmBeta, newdata=testing, type = "response") } StackedPredictions <- data.frame(predictions)  #####Test the number of iterations needed to produce a consistent result dummyM <- data.frame() #dataframe to store mean and median estimates for each iteration IterationsRmse <- data.frame() #dataframe to store rmse of iteration for mean and median for(i in 1:(ncol(StackedPredictions))) {   dummyM[1:nrow(StackedPredictions), "Avg"] <- rowMeans(as.matrix(StackedPredictions[,1:i]))   dummyM[1:nrow(StackedPredictions), "Median"] <- rowMedians(as.matrix(StackedPredictions[,1:i]))      dummyM[1:nrow(StackedPredictions), "AvgError"]<-testing$PercentGravel-dummyM$Avg   dummyM[1:nrow(StackedPredictions), "MedianError"]<-testing$PercentGravel-dummyM$Median      IterationsRmse[1,i] <- sqrt(mean((dummyM$AvgError)^2))   IterationsRmse[2,i] <- sqrt(mean((dummyM$MedianError)^2)) }  it <- data.frame(t(IterationsRmse)) it[3] <- 1:nrow(it) names(it)[1:3] <- c("Mean", "Median", "Iteration")  ggplot(data = it, aes(x = Iteration)) +    geom_line(aes(y = Mean, color = "Mean"), size = 0.5) +    geom_line(aes(y = Median, color = "Median"), size = 0.5) +   labs(title = "Effects of the Number of Bagging Iterations on Root Mean Squared Error", y = "Root Mean Squared Error") +   scale_color_discrete(name="Aggregation Type")  #######All Sand + Gravel Sand100 <- list.files("D:/PeterWhitman_Data/Data/MScResearch/HistEqualizedPhotos/Sand100") Sand <- data.frame(substr(Sand100, 1, nchar(Sand100) - 4)) Sand[, "Estimate"] <- 0 names(Sand)[1:2] <- c("Name", "Estimate")  Gravel100 <- list.files("D:/PeterWhitman_Data/Data/MScResearch/HistEqualizedPhotos/Gravel100") Gravel <- data.frame(substr(Gravel100, 1, nchar(Gravel100) - 4)) Gravel[,"Estimate"] <- 1 names(Gravel)[1:2] <- c("Name", "Estimate")  Percent100 <- rbind(Sand, Gravel)  ####Final Model library(glmmTMB)   69 master <- read.csv("D:/PeterWhitman_Data/Projects/MScResearch/CroppedPhotos/GlCM_Stats/MasterGLCM.csv") master$Brightness <- as.numeric((master$Green + master$Red + master$Blue)/3)  df <- read.csv("TextureStatistics.csv") df$Brightness <- as.numeric((df$Green + df$Red + df$Blue)/3) train <- df %>%   filter(PercentGravel > 0.01 & PercentGravel < 0.99)   AntiMaster <- anti_join(master, Percent100, by = "Name")  AntiPrediction <- data.frame() for(i in 1:1000) {   positions <- sample(nrow(train),size=floor((nrow(train)/4)*3))   training<- train[positions,]      GlmmBeta <- glmmTMB(PercentGravel ~ Entropy + Correlation + Brightness + (1|Day), data = training, family = beta_family(link = "logit"))   AntiPrediction[1:nrow(AntiMaster),i] <- predict(GlmmBeta, newdata=AntiMaster, type = "response") } AntiMaster$Estimate <- rowMedians(as.matrix(AntiPrediction[,1:ncol(AntiPrediction)])) RefinedMaster <- AntiMaster[,c("Name", "Estimate")]  CombinedMaster <- rbind(RefinedMaster, Percent100)  write.csv(CombinedMaster, "GravelEstimates.csv")   Appendix B   This Appendix outlines complementary material for the methods used to analyze and process archived aerial photographs.  B.1 Air Photo Preprocessing Each air photo analyzed in this study included a black or white margin. To ensure that these margins did not confound subsequent analyses, an automated image processing routine was developed in MATLAB. This code segments the scanned air photos to isolate the photo area as a binary image object. The center of this image object should be in the same approximate location as the photo’s principal point. The automated routine then crops the original air photo by a specified pixel distance from the center of the segmented image object. The specified pixel distance for the air photos was 9750, which produced cropped images that were 19501x19501 pixels. For 10 x 10 in contact prints scanned at 2400 dpi these dimensions ensure that the   70 cropped image contains only the photo area while accounting for small amounts of variation in the automated routine. Figure B.1 provides a visual representation of the routine.      The aerial photographs from the 1951 survey contained a reticle representing the location of the principle point. The distance between the principal point identified by the automated routine and the location of the actual principal point was measured within ten photos from the 1951 survey. These measurements are outlined in Table B.1. The average error of the routine within the measured photographs is ~100 pixels with a standard deviation of ~25 pixels.   Table B.1 The distance in pixels and millimeters between the principal point identified by the image processing routine and the actual principal point in ten photos from the 1951 aerial survey. Photo Distance (px) Distance (mm) 1 87.7268488 0.929904597 2 142.3376268 1.508778844 3 104.995238 1.112949523 4 77.07788269 0.817025557 5 98.23441352 1.041284783 6 53.85164807 0.57082747 7 103.8701112 1.101023179 8 92.08691547 0.976121304 9 102.1077862 1.082342534 10 132.2308587 1.401647103 Figure B.1 A visual representation of the automated image processing routine used to crop each archived photograph. The original image (left) is segmented to identify the photo area and by extension the principle points (middle). The original image is cropped by a specified distance from the identified principal point (right).   71 The MATLAB script developed to crop the photographs from each aerial survey:  imageDir = %Input folder directory. Folder name should end with %the year in which the survey took place (e.g., AirPhotos1951) TestDir =  %Output folder directory filePattern = fullfile(imageDir); ImageFiles = dir(filePattern); series = imageDir((length(imageDir)-3):length(imageDir));  for i = 3:length(ImageFiles) baseFileName = ImageFiles(i).name; fullFileName = fullfile(imageDir, baseFileName); I = imread(fullFileName);  %The photographic paper used for the 1951 photos is very different than the %photographic paper used for the photos taken in 2008, 1999, and 1984. %Therefore, the 1951 photos must be processed with different parameters t = strcmp(series, '1951'); if t == 1     BW1 = imbinarize(I, 0.9); else     BW1 = imbinarize(I, 0.3); end  %Segementation routine to isolate the margins from the photo area BW2 = imcomplement(BW1); BW3 = bwareafilt(BW2, 1); BW4 = imcomplement(BW3); BW5 = bwareafilt(BW4, 1);      if t == 1     BW6 = imcomplement(BW5);     BW_stat = regionprops(BW6,'Centroid'); else     BW_stat = regionprops(BW5,'Centroid'); end  %These dimensions should be good for all 10"x10" photos scanned at 2400 dpi width = 19500; height = 19500;  %Location of the pricipal point (nadir) upper_x = abs((BW_stat.Centroid(:,1) - (width/2))); upper_y = abs((BW_stat.Centroid(:,2) - (height/2)));  croppedI = imcrop(I, [upper_x, upper_y, 19500, 19500]); TestSave = fullfile(TestDir, baseFileName); imwrite(croppedI, TestSave); end         72 B.2 Orthomosaic Workflow  Figure B.2 A flowchart of the orthomosaic workflow carried out in Agisoft Photoscan                            73  Table B.2 The input parameters used to construct the archived and contemporary orthomosaics in Agisoft Photoscan   Archived Contemporary General   Markers None 9 Coordinate system Local Coordinates (m) WGS 84 / UTM zone 10N Rotation angles Yaw, Pitch, Roll Yaw, Pitch, Roll Point Cloud   Accuracy High Highest Generic preselection Yes Yes Key point limit 400,000 40,000 Tie point limit 40,000 4,000 Adaptive camera model fitting Yes Yes Parameters f, b1, b2, k1, k2, p1, p2 f, cx, cy, k1, k2 Depth Maps   Quality Medium Medium Filtering mode Aggressive Aggressive Dense Point Cloud   Quality Medium Medium Depth filtering Aggressive Aggressive 3D Model   Surface type Height field Height field Source data Dense Dense Interpolation Extrapolated Enabled Quality Medium Medium Depth filtering Aggressive Aggressive Orthomosaic   Coordinate system Local Coordinates (m) WGS 84 / UTM zone 10N Colors 1 bands, uint8 4 bands, uint8 Blending mode Mosaic Mosaic Surface Mesh DEM Enable hole filling Yes Yes    74 B.3 Merge and Scale Parameters Figure B.3 The overall accuracy produced by each combination of merge and scale parameters at intervals of ten when used to classify the contemporary orthomosaic after it was matched to the 1951 aerial survey. Random trees classifier (top) and support vector machine (bottom).   75   Figure B.4 The overall accuracy produced by each combination of merge and scale parameters at intervals of ten when used to classify the contemporary orthomosaic after it was matched to the 1984 aerial survey. Random trees classifier (top) and support vector machine (bottom).   76   Figure B.5 The overall accuracy produced by each combination of merge and scale parameters at intervals of ten when used to classify the contemporary orthomosaic after it was matched to the 1999 aerial survey. Random trees classifier (top) and support vector machine (bottom).   77  Figure B.6 The overall accuracy produced by each combination of merge and scale parameters at intervals of ten when used to classify the contemporary orthomosaic after it was matched to the 2008 aerial survey. Random trees classifier (top) and support vector machine (bottom). 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0380528/manifest

Comment

Related Items