UNIVERSIDAD COMPLUTENSE DE MADRID FACULTAD DE GEOGRAFÍA E HISTORIA MÁSTER UNIVERSITARIO EN TECNOLOGÍAS DE LA INFORMACIÓN GEOGRÁFICA TRABAJO FIN DE MÁSTER CURSO 2019-2020 Empleo de técnicas de fotogrametría para obtener modelos digitales de superficie a partir de fotografías aéreas históricas y detectar cambios geomorfológicos: aplicación al valle del río Guadalfeo (provincia de Granada). Use of photogrammetry techniques to obtain digital surface models from historical aerial photographs and detect geomorphological changes: application to the valley of the Gualdafeo river (Granada province). David Nievas Valverde CONVOCATORIA: SEPTIEMBRE TUTOR: Luis Miguel Tanarro Garcia Departamento de Geografía. Universidad Complutense de Madrid 2 RESUMEN El presente Trabajo Fin de Máster (TFM) toma como área de estudio una sección del tramo medio del Río Guadalfeo, al sur de la provincia de Granada y muestra la metodología llevada a cabo mediante el uso de las TIGs, y en especial la aplicación de técnicas fotogramétricas “Structure from Motion-MultiView Stereo” (SfM-MVS), con la finalidad de obtener modelos 3D del área de estudio a partir de fotogramas aéreos históricos (vuelo Americano Serie B, de 1956-57), así como a partir de fotogramas más recientes pertenecientes a los vuelos del PNOA de los años 2006 y 2016 para la obtención de productos cartográficos secundarios. Estos productos secundarios consisten en una ortofoto, un modelo digital de superficies (MDS) y una nube de puntos densa para cada año estudiado. Posteriormente, a partir los productos obtenidos, se ha procedido en primer lugar a realizar una cartografía de la llanura de inundación y sus principales elementos geomorfológicos, incluyendo una serie de deslizamientos presentes en las laderas del valle. En segundo lugar, mediante el software de código abierto “Geomorphic Change Detection” (GCD), se ha procedido a generar dos modelos digitales de diferencias de elevación (DoD) mediante la resta de los MDS de los años 1956-2006 y 2006-2016. Por último, dichos productos han sido utilizados para su correspondiente evaluación y comparación, permitiendo analizar los cambios producidos en la llanura de inundación a lo largo del intervalo temporal. Por un lado, la observación de las ortofotos ha permitido mostrar la variabilidad geomorfológica del canal o lecho del río; y, por otro lado, los resultados de los DoD han permitido reflejar los tramos de la llanura de inundación donde han ocurrido cambios asociados a procesos de sedimentación o erosión o al crecimiento de la vegetación. Palabras clave: Structure from Motion-MultiView Stereo (SfM-MVS); fotogrametría digital, DoD (Modelos Digitales de Diferencias de Elevación); geomorfología fluvial, llanura de inundación, río Guadalfeo. ABSTRACT The present Master Thesis takes as a study area a section of the middle stretch of the Guadalfeo River, south of the province of Granada and shows the methodology carried out through the use of TIGs, and especially the application of photogrammetric techniques Structure from Motion-MultiView Stereo (SfM-MVS), in order to obtain 3D models of the area from historical aerial photographs (American Flight Series B, 1956- 57) as well as from more recent photographs from PNOA flights in 2006 and 2016 for the obtaining of secondary cartographic products. These secondary products consist of an orthophoto, a digital surface model (DSM) and a dense point cloud for each year studied. Subsequently, based on the products obtained, a cartography of the floodplain and its main geomorphological features was carried out, including a series of landslides present in the slopes of the valley. Second, using the Geomorphic Change Detection “GCD” open-source software, two digital elevation difference models (DoD) have been generated by subtracting the MDS for the years 1956-2006 and 2006-2016. 3 Finally, these products have been used for their evaluation and comparison, making it possible to analyse the changes produced in the floodplain throughout the time interval. On the one hand, the observation of the orthophotos has shown the geomorphological variability of the channel or riverbed; and, on the other hand, the DoD results have made it possible to reflect the stretches of the floodplain where changes associated with sedimentation or erosion processes or the growth of vegetation have occurred. Keywords: SfM-MVS; digital photogrammetry, DoD (DEMs of Difference); fluvial geomorphology, floodplain, Guadalfeo river. 4 ÍNDICE 1. INTRODUCCIÓN .................................................................................................... 5 1.1. Importancia de las Tecnologías de la los Sistemas de Información Geográfica (TIG) en el análisis geomorfológico. ............................................................................ 5 1.2. La técnica fotogramétrica “Structure from Motion” (SfM) y la obtención de Modelos Digitales de Superficies (MDS) y otros productos (ortofotos, nubes de puntos) .......................................................................................................................... 6 1.3. Los Modelos Digitales de Superficies (MDS) y los modelos de diferencias de elevación o superficies (DoD). ..................................................................................... 7 2. OBJETIVOS............................................................................................................ 10 3. ÁREA DE ESTUDIO: contexto geográfico y geológico. ....................................... 11 4. MATERIALES Y METODOLOGÍA ..................................................................... 14 4.1 Obtención de materiales ........................................................................................ 14 4.2. Procesamiento de las fotografías aéreas usando fotogrametría digital automática o Structure from Motion. ............................................................................................... 15 4.2.1. Obtención de puntos de control en los fotogramas de cada vuelo. ............... 16 4.3. Aerotriangulación o alineación de las fotografías aéreas o fotogramas. ............. 22 4.3.1. Control de calidad del proceso de aérotriangulación en ContextCapture. .... 24 4.4. Producción del modelo 3D y de otros productos cartográficos. .......................... 25 3.4.1. Productos secundarios del procesamiento fotogramétrico. ........................... 27 4.5. Reconocimiento y digitalización de los principales elementos geomorfológicos sobre las ortofotos. ...................................................................................................... 29 4.6. Método para la elaboración de un Modelo Digital de diferencias de Superficie con el software Geomorphic Change Detection (GCD). ............................................ 30 5. RESULTADOS ....................................................................................................... 37 5.1. Análisis geomorfológico de la evolución del río Guadalfeo a partir de la interpretación de las ortofotos. .................................................................................... 37 5.2. Análisis de los cambios de elevación o alturas a partir de DoDs. ....................... 41 5.2.1 Análisis de los cambios de elevación o alturas en la llanura de inundación para los años 1956 y 2006. ...................................................................................... 41 5.2.2 Análisis de los cambios de elevación o alturas en la llanura de inundación para los años 2006 y 2016. ...................................................................................... 44 5.2.3 Análisis de los cambios de elevación o alturas en los deslizamientos de la ladera del valle del río Guadalfeo para los años 1956 y 2006. ................................ 47 5.2.4. Análisis de los cambios de elevación o alturas en los deslizamientos de la ladera del valle del río Guadalfeo para los años 2006 - 2016. ................................ 50 6. CONCLUSIONES .................................................................................................. 53 7. BIBLIOGRAFÍA ..................................................................................................... 55 5 1. INTRODUCCIÓN 1.1. Importancia de las Tecnologías de la los Sistemas de Información Geográfica (TIG) en el análisis geomorfológico. Los “Sistemas de Información Geográfica” (SIG), dentro de las TIG, son considerados como una herramienta fundamental para lograr la integración de diferentes tipos de datos de diversos orígenes y fuentes debido al hecho de que estos comparten un mismo ambiente o entorno, es decir presentan un marco geográfico y un “universo digital”, que generalmente se expresa como la pantalla de una computadora. Apoyándose en el principio de superposición de capas y en el manejo de datos geográficos o espaciales, los SIG permiten almacenar, recuperar, manipular, analizar e identificar relaciones espaciales en forma de mapas o datos descriptivos. (Asato et al., 1996). En las disciplinas vinculadas con las Ciencias de la Tierra, como la geomorfología, los SIG se han utilizado para desarrollar estudios diversos. Cabe mencionar, por ejemplo, la serie de publicaciones del Instituto Geológico y Minero de España sobre mapas de susceptibilidad a los movimientos de ladera con técnicas SIG (Ayala-Carcedo y Corominas, 2002), los sistemas de información geográfica en la gestión de los riesgos geológicos y el medioambiente (Laín-Huerta, 2002) o en la elaboración de mapas de peligrosidad de avenidas e inundaciones (Díez-Herrero et al., 2006.). En estos libros se presentan y se resuelven problemas de índole geomorfológico mediante la utilización de herramientas SIG. Además, los SIG permiten contar con resultados representados en forma gráfica, consiguiéndose una rápida evaluación de los proyectos y facilitando su compresión por parte de personas ajenas al ambiente geomorfológico (Asato et al., 1996). Desde esta perspectiva, este Trabajo Fin de Máster (TFM) se presenta con la motivación de plantear un problema geomorfológico asociado con la detección de cambios en la superficie topográfica entre varias fechas, es decir cambios temporales, que se manifiestan en erosión y/o acumulación de materiales. Los SIG, y otras tecnologías como la fotogrametría, se han empleado con frecuencia para estudiar los procesos de erosión-sedimentación, principalmente, en áreas acarcavadas (Nadal et al., 2015), pero también en lechos fluviales (Cook, 2017; Llena et al., 2018; Nourbakhshbeidokhti et al., 2019). No obstante, el análisis para detectar cambios geomorfológicos temporales, en especial en periodos cortos, precisa disponer de modelos digitales de elevaciones (MDE) y ortofotos de alta resolución, con tamaño de celda centimétrica o milimétrica. La aparición reciente de herramientas geoespaciales (LíDAR, TLS, SfM) permite la obtención de modelos topográficos tridimensionales de alta resolución. Asimismo, estas herramientas ofrecen la posibilidad, también, de crear ortofotos de alta calidad y resolución. Entre las técnicas citadas, la técnica fotogramétrica “Structure from Motion (SfM)” se puede aplicar a fotografías aéreas históricas, cuyo resultado puede ser la obtención de Modelos Digitales de Superficies (MDS) y ortofotos de muy alta resolución (Gomez, 2012; Lane et al., 2010; Llena et al., 2018). A partir de estas consideraciones, este TFM persigue, como más adelante se explicará con más detalle, obtener MDS y ortofotos de fotografías aéreas históricas (vuelo 6 americano de 1956/57) y de fotografías digitales del Plan Nacional de Ortofotografía Aérea (PNOA) de los años 2006 y 2016 mediante técnicas fotogramétricas, en un tramo del valle del río Guadalfeo (Granada). Así, los MDS y ortofotos obtenidas se podrán comparar para identificar y cuantificar los cambios geomorfológicos. Previamente a explicar los objetivos y la metodología se expone, brevemente, los fundamentos teóricos de la técnica fotogramétrica “Structure from Motion (SfM) utilizada para obtener MDS, así como el método más usual para proceder a su comparación. 1.2. La técnica fotogramétrica “Structure from Motion” (SfM) y la obtención de Modelos Digitales de Superficies (MDS) y otros productos (ortofotos, nubes de puntos) En los últimos años, el desarrollo de las tecnologías fotogramétricas de foto- reconstrucción 3D, tales como “Structure from Motion” junto con la aplicación de algoritmos “Multiview-Stereo (MVS)” han mejorado la calidad y resolución de los datos 3D obtenidos del terreno. Debido a que en la gran mayoría de estudios la implementación de los algoritmos MVS está presente, lo más correcto sería hacer referencia a este conjunto de procesos como “SfM-MVS” (Micheletti et al., 2015; Smith et al, 2015). La técnica “SfM-MVS” se puede definir como un método fotogramétrico digital y automatizado de alta resolución y bajo coste capaz de componer una estructura 3D a partir de la superposición o solapamiento de imágenes obtenidas desde distintos puntos de vista (Tomás et al., 2016; Muñoz-Narciso et al., 2014). Es decir, gracias a esta técnica podemos construir estructuras tridimensionales a partir de imágenes bidimensionales por medio de la extracción de puntos comunes entre los distintos fotogramas obtenidos. Esta identificación y extracción de puntos comunes es llevada a cabo de forma general mediante algoritmos “SIFT” o “Scale Invariant Feature Transform”. A partir de los puntos comunes extraídos se orientan y localizan las cámaras y, posteriormente, se obtiene una nube de puntos de baja densidad en tres dimensiones (o malla o modelo 3D) en un sistema de coordenadas local. Posteriormente, dicha nube de puntos es intensificada mediante el uso de los mencionados algoritmos “MVS” con el fin de generar una nube de puntos densa a partir de la cual es posible elaborar un mosaico orto-referenciado de imágenes, modelos digitales superficies (MDS) de alta resolución, contornos topográficos y modelos 3D (Fernández-Lozano y Gutiérrez-Alonso, 2016; Micheletti et al., 2015). En contraposición con la fotogrametría tradicional, “SfM-MVS” no requiere datos de orientación ni ubicación 3D de las cámaras para la reconstrucción de la escena ya que las posiciones de cámara y la orientación se resuelven automáticamente y sin la necesidad de establecer una red de puntos de control de coordenadas 3D conocidas (Tomás et al., 2016; Muñoz-Narciso et al., 2014). Sin embargo, el uso de puntos de control nos permitirá no sólo referir nuestro modelo 3D a un sistema de coordenadas geográfico conocido, sino que también será de utilidad a la hora minimizar los errores de reproyección y georreferenciación, permitiéndonos, por tanto, obtener una mejor 7 estimación de los parámetros intrínsecos de la cámara y una mejor reconstrucción de la escena (Muñoz-Narciso et al., 2014; Smith et al., 2015) En fecha reciente son numerosos los estudios que han hecho uso de las mencionadas técnicas a fin de generar productos relacionados con la topografía (Fonstad et al., 2013; Gomez, 2012; James y Robson, 2012). Además, dichos estudios han mostrado como “SfM-MVS” garantiza un bajo coste a la hora de obtener datos topográficos, lo cual es de gran interés especialmente en aquellas investigaciones dónde los cambios en el terreno se producen con una considerable rapidez (Muñoz-Narciso et al., 2014). Adicionalmente y de forma contraria a las técnicas tradicionales fotogramétricas, donde se requiere un equipo costoso y un usuario especializado y experimentado en la materia, el desarrollo de los métodos “SfM-MVS” ofrece la oportunidad de adquirir a muy bajo costo datos tridimensionales con una menor supervisión y experiencia por parte del usuario (Westoby et al., 2012; Micheletti et al., 2015). 1.3. Los Modelos Digitales de Superficies (MDS) y los modelos de diferencias de elevación o superficies (DoD). Cuantificar los cambios volumétricos es una de las prioridades en muchos de los estudios en geomorfología, enfocados a estudiar los cambios sobre el terreno o la superficie topográfica. Gracias a que en los últimos años tanto las técnicas de procesamiento de datos topográficos como las tecnologías geoespaciales han experimentado un gran desarrollo, ha sido posible en el ámbito de la geomorfología llevar a cabo campañas de monitoreo capaces de obtener modelos digitales precisos y de alta resolución, acordes con las tasas de evolución o cambio de la forma del relieve. Además, las nuevas tecnologías de análisis y procesado de imágenes nos ofrecen la posibilidad de trabajar con imágenes históricas pudiendo obtener modelos digitales a partir de estas también. Este hecho brinda la oportunidad de extender las escalas de tiempo de las investigaciones, así como realizar inferencias sobre los procesos que pudieron ocasionar determinados cambios sobre el terreno, basadas en los datos obtenidos durante los periodos de monitoreo (Williams, 2012). Entre los productos generados a partir de la técnica “Structure from Motion” cabe destacar la obtención de Modelos Digitales de Superficies (MDS), que presenta algunas diferencias con los denominados Modelos Digital de Elevaciones (MDE). La diferencia entre ambos tipos de modelos radica fundamentalmente en que mientras el MDS representa el relieve incluyendo todo tipo de elementos presentes como edificaciones, vegetación u otros objetos, el MDE representa únicamente el suelo desnudo. Como desventaja, dicha diferencia puede llegar a suponer a veces un problema ya que con un MDS pueden generarse problemas tales como vacíos en la interface tierra-agua, sombreado e inversión de relieve, así como problemas en el modelado hidrológico (Correa-Muñoz, 2012). En todo caso, la disponibilidad de MDE (o MDS) es una variable fundamental en cualquier estudio geomorfológico, ya que proporcionan una estructura numérica de datos que representa la distribución espacial de la altitud de la superficie del terreno (Felicísimo, 1994). En otras palabras, gracias a dicha técnica podremos obtener una imagen ráster en la cual cada píxel estará representado por un valor de altitud. 8 A raíz de la presencia de estos modelos digitales de elevaciones y/o superficies, hoy en día es también posible generar mediante determinadas herramientas los llamados “DEMs of Difference (DoDs)”. El acrónimo DoD hace referencia a un modelo digital de diferencias de elevación que es el resultado de restar los valores de elevación de un MDE más antiguo con los de uno más reciente, celda por celda (http://gcd.riverscapes.xyz/Concepts/dod.html). Es decir, de esta manera obtendremos un nuevo MDE a partir de la resta de otros dos en el cual los valores positivos representarán aquellas zonas en las que ha habido sedimentación y valores negativos en las zonas donde se ha producido erosión dentro del intervalo de tiempo comprendido entre la toma de los datos que se usaron para generar el modelo más antiguo y el más reciente (Figura 1). Mientras que estos productos son de gran utilidad para el monitoreo y la interpretación geomorfológica, también es cierto que los errores en los datos y modelos pueden propiciar interpretaciones erróneas, por lo que es de gran importancia realizar un análisis previo de los posibles errores, para tenerlos en cuenta a la hora de generar el DoD, e interpretar los resultados correctamente. Es por ello que a su vez se están desarrollando nuevos métodos con los que es posible estimar de una forma más fiable las incertidumbres en los modelos digitales con el fin de obtener una mejor estimación de los cambios geomorfológicos (Wheaton et al., 2010) Figura 1. Ejemplo de flujo de trabajo para obtener un “DoD (DEM of Difference). Fuente: “http://gcd.riverscapes.xyz/Concepts/challenges.html” 9 Sin duda, el desarrollo y avance de las nuevas tecnologías (LiDAR, SfM-MVS) combinado con la obtención de fotografías por medios más sencillos y de menor coste (drones, cámaras comerciales, etc.) permite obtener MDS o MDE en diferentes fechas que, posteriormente pueden ser comparados. En este sentido, los DoDs se están aplicando en geomorfología fluvial en aspectos tales como para calcular las tasas de transporte de sedimentos, para la interpretación de procesos fluviales, como zapamiento lateral, llenado, migración y avulsión de canales, entre muchos otros (Muñoz–Narciso et al., 2014; Williams, 2012; Cook, 2017); o incluso para detectar los cambios espaciales en intervalos de tiempo más prolongados a partir del procesamiento de fotografías aéreas historias con técnicas fotogramétricas, como SfM-MVS (Lane et al., 2010; Llena et al., 2018). Por supuesto, en el caso de este TFM, esta técnica supone una excelente aproximación a la hora de ayudarnos a determinar cambios en la dinámica geomorfológica del tramo fluvial estudiado, o más bien en qué zonas del mismo se ha producido sedimentación y en cuales erosión a lo largo del tiempo. Todo ello mediante la resta de los distintos MDS que serán obtenidos para distintas fechas. 10 2. OBJETIVOS El objetivo principal de este TFM es valorar la técnica fotogramétrica “Structure from Motion” (SfM), aplicada a fotografías aéreas, para estudiar la dinámica fluvial de un tramo del río Guadalfeo (Granada). Para ello, se han seleccionado tres vuelos: 1956/57 (vuelo americano, serie B), 2006 y 2016 (vuelos del PNOA), que cubren un periodo de unos 60 años. Las fotografías se han procesado en un programa de fotogrametría comercial (Bentley ContextCapture), lo que ha permitido obtener el Modelo Digital de Superficies y la ortofoto de cada vuelo. Posteriormente, los MDS se han comparado para generar Modelos Digitales de Diferencias de Elevación (que como se ha explicado se citan en la bibliografía con el acrónimo DoD, en inglés DEM of Differencing) con el fin de detectar los cambios geomorfológicos en términos de erosión y sedimentación, que han ocurrido entre estas fechas en la llanura de inundación del río, y en algunos sectores de ladera afectados por deslizamientos. Para cumplir este objetivo general, se han cubierto los siguientes objetivos secundarios: - Identificación de puntos de control (PC) en los fotogramas de cada vuelo para obtener sus coordenadas “x,y” sobre la ortofoto y “z” sobre un MDE, creado con el vuelo LíDAR (1ª Cobertura). - Georreferenciación y procesamiento (aerotriangulación) de las fotografías aéreas de 2016, 2006 y 1956 en el programa de fotogrametría, y determinación de los errores de la georreferenciación. - Producción del modelo 3D de cada vuelo y generación de los productos cartográficos derivados (MDS, ortofotos y nubes de puntos densas) como productos secundarios de las foto-reconstrucciones 3D. - Comparación de los MDS generados mediante la obtención de Modelos Digitales de Diferencias de Elevación, que ofrecerán valores negativos (erosión) o positivos (sedimentación) (Williams, 2012). Para ello se ha utilizado el programa específico de código abierto GCD (acrónimo en ingles de Geomorphic Change Detection, http://gcd.riverscapes.xyz/), diseñado para sustraer MDS o MDE. El resultado permitirá comparar los resultados de erosión-sedimentación y volumen obtenidos en las fechas estudiadas. - Observación y digitalización sobre las ortofotos de las principales formas que se reconocen en la llanura de inundación del río con el fin de elaborar mapas geomorfológicos sencillos. Su observación permitirá mostrar la dinámica de algunos elementos geomorfológicos, como son las variaciones del lecho fluvial. 11 3. ÁREA DE ESTUDIO: contexto geográfico y geológico. El área de estudio, como se ha dicho, se corresponde con un tramo del valle del río Guadalfeo, al sur de la provincia de Granada (Figura 2). Este río, cuya cuenca hidrográfica ocupa una amplia extensión (1.295 km2), tiene su nacimiento en las estribaciones de la vertiente meridional de Sierra Nevada, situándose su cabecera a cotas próximas a los 3000 metros (Peñón del Puerto, 2.754 m). La cuenca hidrográfica de este río, que se asienta, en gran parte, sobre la comarca de Las Alpujarras, con altitudes que superan los 2.000 m.s.n.m, es drenada por importantes ríos de montaña que descienden desde las cumbres (Mulhacén a 3479 m, Veleta a 3398) y altas laderas de Sierra Nevada, como, por ejemplo, los ríos Lanjarón, Poqueira, Trevelez. El río Guadalfeo, que recorre aproximadamente 71 km hasta desembocar en las proximidades del núcleo de población de Salobreña, presenta un régimen irregular de carácter nivo- pluvial, pues su caudal varía considerablemente durante el año de unas estaciones a otras. En concreto, el tramo seleccionado se encuentra localizado en un sector de la cuenca media, concretamente entre la localidad de Torvizcón y aguas arriba de la presa de Rules, que fue inaugurada en el año 2004. Desde un punto de vista geológico, la consulta de la hoja y memoria 1042 (Lanjarón) del Mapa Geológico de España a escala 1:50.000, las laderas y fondo del valle están Figura 2. Mapa de la cuenca hidrográfica del río Guadalfeo. El rectángulo rojo indica el tramo fluvial seleccionado. 12 constituidas por materiales litológicos correspondientes al complejo Alpujárride, además de depósitos neógenos y cuaternarios de carácter postorogénico. Dentro del área estudiada, los materiales del complejo Alpujárride, que se sitúa dentro de las zonas internas de la Cordillera Bética entre los complejos Nevado-Filábride y Maláguide, se distribuyen en tres unidades tectónicas o mantos (IGME, 1979; Jiménez Perálvarez, 2005): • Manto del Alcazar. Dentro del área de estudio, las formaciones que afloran están constituidas por filitas y cuarcitas, y localmente yesos y calcoesquistos, de edad Triásico Inferior-Pérmico. • Manto de Cástaras. Las formaciones principalmente representadas en el tramo del valle estudiado están formadas por calizas recristalizadas y dolomías de Triásico Medio-Superior. • Manto de Lújar. Los materiales que afloran, dentro las formaciones de este manto, en el área estudiada se corresponden con calizas y dolomías con intercalaciones locales de margocalizas, yeso y arcillitas, también del Triásico Medio-Superior. Los materiales néogenos, que aparecen ampliamente en un sector del margen derecho del valle, descansando sobre los materiales de los mantos alpujárrides, están constituidos por conglomerados heterométricos con intercalaciones arenosas del Plioceno. Finalmente, los depósitos cuaternarios se corresponden con aluviones actuales (cauce o canal) o recientes (terrazas bajas y/o lecho mayor) formando la llanura de inundación del río Guadalfeo y de los principales arroyos. Además, junto a estas formaciones, afloran conglomerados cementados asociados a depósitos de ladera o abanicos aluviales. Sobre estos materiales litológicos, el río Guadalfeo se ha encajado profundamente y ha modelado laderas muy escarpadas donde los procesos geomorfológicos son muy activos (Galve et al., 2017). Así, estas laderas de fuerte pendiente están afectadas por procesos de inestabilidad, dando lugar a frecuentes movimientos en masa (desprendimientos, deslizamientos, debris flow) que afectan principalmente a las laderas de la margen norte del río Guadalfeo (Tomás Fernández et al., 2017; Palenzuela et al., 2016; Fernández Chacón et al., 2015; Jiménez Perálvarez et al., 2010). Asimismo, la evolución del lecho fluvial del río Guadalfeo bajo unas condiciones medioambientales semiáridas ha propiciado, en su tramo medio, el desarrollo de un lecho fluvial de tipo “rambla”, que se caracteriza por un fondo plano en el cual alternan corrientes efímeras con episodios de inundación. Éstos están controlados por tormentas intensas o por la fusión rápida de la nieve de Sierra Nevada (Tomás Fernández et al., 2017). De este modo, el río Guadalfeo, entre aguas arriba de la presa de Rules y el paraje de la Sierra de Jubiley, presenta un modelado característico donde se reconocen elementos y rasgos geomorfológicos de un curso de agua de tipo “rambla o uadi” (Coque, 1987). Los episodios de inundación son frecuentes y provocan la alteración y modificación de la llanura de inundación debido a la migración del canal o cauce hacia las orillas (Figura 3b) (Fernández-Chacón et al., 2015). Teniendo en cuenta estos aspectos geomorfológicos el propósito de este TFM, como se ha dicho, se centra, principalmente, en estudiar la dinámica de esta llanura de 13 inundación, y secundariamente la de los deslizamientos, a través del empleo fotografías aéreas de varias fechas. Figura 3a. Fotografía de un tramo del Río Guadalfeo. Fuente: “https://es.wikiloc.com/rutas-senderismo/pitres- orgivapasando-por-finca-el-valero-11441089/photo-6992031”. Figura 3b. Fotografía de un tramo del Río Guadalfeo tomada durante una inundación el día 24 de diciembre de 2009. Fuente: “https://granadanatural.com/ficha_paisajes.php?cod=263”. https://es.wikiloc.com/rutas-senderismo/pitres-orgivapasando-por-finca-el-valero-11441089/photo-6992031 https://es.wikiloc.com/rutas-senderismo/pitres-orgivapasando-por-finca-el-valero-11441089/photo-6992031 https://granadanatural.com/ficha_paisajes.php?cod=263 14 4. MATERIALES Y METODOLOGÍA 4.1 Obtención de materiales Para la realización de este estudio ha sido necesario recopilar los materiales que se sintetizan en la Tabla 1. Los materiales de partida fundamentales han sido las fotografías aéreas de las tres fechas que se han seleccionado para este estudio. Así, los fotogramas, correspondientes al vuelo americano de 1956-57 (serie B), se compraron en el Centro Cartográfico y Fotográfico (CECAF) - Ejército del Aire. Para cubrir el área de estudio fue preciso adquirir 9 fotogramas, los cuales son escaneados en el CECAF a 21 micras (1209 ppp). Por su parte, los fotogramas más recientes (vuelos del PNOA de 2006 y 2016) se descargaron del portal web oficial del Centro Nacional de Información Geográfica (CNIG) del Instituto Geográfico Nacional (IGN). Se precisaron 16 y 17 fotogramas respectivamente, y a diferencia del vuelo anterior, están en formato digital, y ambos tienen una resolución de 0,45 m También se descargaron del portal del CNIG-IGN, los datos del vuelo LIDAR (1ª Cobertura) con el fin, como luego se explicará, de crear el modelo digital de elevaciones como base para obtener la coordenada Z de los puntos de control para georreferenciar los fotogramas de los vuelos; y, con mismo fin se usó la ortofoto de 2016, pero en este caso para obtener las coordenadas X e Y. Material y año Proyecto Nº Fotogramas / Formato Escala o Resolución Fuente Fotogramas 1956 Vuelo Americano Serie B (1956- 1957) 12 / *jpg 1:33.000 (fotogramas analógicos escaneados a 21 micras) Centro Cartográfico y Fotográfico (CECAF) - Ejército del Aire Fotogramas 2006 Plan Nacional de Ortofotografía Aérea PNOA (PNOA) 16 / *ecw 0,45 m Centro Nacional de Información Geográfica (CNIG) - Instituto Geográfico Nacional (IGN) Fotogramas 2016 17 / *ecw 0,45 m LIDAR 1ª Cobertura (2008-2015) *LAZ 2x2 km Ortofoto 2016 *ecw 0,50 m Además, en este Trabajo Fin de Máster se han utilizado otros materiales complementarios, que se han empleado como referencia para completar diferentes apartados de la Memoria. Entre estos cabe destacar, la serie de mapas topográficos a escala 1:50.000 y 1:25.000, y la hoja 1042 (Lanjarón) del mapa geológico a escala 1:50.000. Tabla 1. Material de partida. 15 Material Formato Sistema Geodésico de Referencia Escala Fuente Mapa Topográfico Nacional, hoja 1042 MTN50 ráster georreferenciado (*tif) UTM ETRS89- 30N 1:50.000 Centro Nacional de Información Geográfica (CNIG) -Instituto Geográfico Nacional (IGN) Mapa Topográfico Nacional, hoja 1042, C-III y IV MTN25 ráster georreferenciado UTM ETRS89- 30N 1.25.000 Mapa Geológico de España, hoja 1042 IGME-50 mapa escaneado (*jpg) y memoria UTM ED50 – Elipsoide Internacional 1:50.000 Instituto Geológico y Minero de España (IGME) A partir de estos materiales, la aplicación de la técnica fotogramétrica y el método de comparación de los DSM y ortofotos, propuestos en este TFM para cumplir con los objetivos, cubrió las siguientes etapas: 4.2. Procesamiento de las fotografías aéreas usando fotogrametría digital automática o Structure from Motion. En este TFM se ha utilizado el software ContextCapture (Bentley), que utiliza un método de restitución basado en algoritmos de tipo “Structure from Motion” y “Multi- View Stereo) (SfM-MVS). El procesamiento y la reconstrucción tridimensional de las fotografías y la producción de modelos 3D y otros productos derivados con el software ContextCapture requiere una serie de pasos consecutivos (Figura 4): añadir las fotografías, especificar las propiedades de la cámara (distancia focal y tamaño del sensor), proceder a alinear o aerotriangular las fotografías, indicar los ajustes de la reconstrucción o precisión geométrica, y, finalmente, obtener los productos cartográficos (modelo o malla 3D, nubes de puntos tridimensionales, MDS, ortofotos) en los formatos, sistema de referencia espacial y resolución espacial que se desee. Aunque previamente antes de realizar la alineación o aerotriangulación de los fotogramas de cada vuelo es posible añadir puntos de control, de modo que los productos generados queden georreferenciados. Tabla 2. Material complementario. Figura 4. Pasos a seguir para la producción de un modelo 3D mediante el software ContextCapture. 16 4.2.1. Obtención de puntos de control en los fotogramas de cada vuelo. Como se acaba de indicar, para el tratamiento y procesamiento de los fotogramas es necesario encontrar o localizar previamente puntos de control con coordenadas XYZ con el fin de que el Modelo 3D resultante de las tres fechas quede correctamente georreferenciado. Estos puntos de control (PC) son puntos que se han buscado de forma manual en los fotogramas y que posteriormente se han localizado también en un mapa de base georreferenciado, sobre el cual dichos puntos de control (PC) se han digitalizado. Ese mapa base ha sido, en nuestro caso, la ortofoto de 2016 del PNOA. Además, en este proceso de identificación de PC, es necesario tener en cuenta que el software con el cual se ha realizado las foto-reconstrucciones 3D de cada vuelo (ContextCapture) requiere que cada punto de control aparezca como mínimo en tres fotogramas distintos. Además de esto, es conveniente que los PCs estén bien distribuidos para que la calidad de la georreferenciación de los modelos sea más precisa. Teniendo en cuenta estas consideraciones el proceso para localizar y obtener puntos de control en los fotogramas ha sido el siguiente: 1. Cargamos los fotogramas del vuelo del PNOA, que tienen una georreferenciación “grosera”, en el software ArcMap junto con la ortofoto de 2016 y establecemos en las propiedades del marco de datos el sistema de coordenadas con el que vamos a trabajar. En este caso: ETRS 1989 UTM Zone 30N. 2. Establecemos unas guías para definir aproximadamente las franjas en las que se solapan tres fotogramas a la vez, con el fin de obtener puntos de control en dichas franjas de forma más sencilla (Figura 5). En el caso de los fotogramas de 1957, éstos no están georreferenciados, por lo que construiremos manualmente un mosaico de los mismos usando el software “Adobe Photoshop CS6” (Figura 6). Posteriormente georreferenciaremos el mosaico con la herramienta “Georreferencing” de ArcMap de forma aproximada con el único fin de tener una guía a la hora de ubicar y localizar los PC comunes en cada fotograma. Figura 5. Uso de guías para definir las franjas en las cuales se solapan tres fotogramas a la vez 17 Cabe destacar que el solapamiento recomendado en las especificaciones del Manual de Bentley Context Capture (FAQS de Bentley) de los fotogramas para obtener una mejor precisión en los resultados es de un 70%, mientras que el mínimo requerido es de un 50%. En nuestro caso, mediante unas mediciones de los fotogramas en ArcGIS y aplicando una regla de tres podemos determinar que tanto los fotogramas de 2006 como los de 2016 presentan un solapamiento entorno al 60% (Figura 7 izquierda). En cuanto a los fotogramas de 1956, con ellos no podemos realizar la misma medición, ya que no se encuentran geo-referenciados ni orto-rectificados, pero haciendo una comparación visual usando Photoshop puede observarse que el solapamiento debe encontrarse también cercano al 60% (Figura 7 derecha). Figura 6. Mosaico de los fotogramas del vuelo de 1956-57 realizado mediante Adobe Photoshop CS6, con el fin de identificar la franja de solapamiento entre tres fotogramas. Figura 7. Las zonas en rojo delimitan el área de solapamiento entre dos fotogramas. La imagen izquierda corresponde a dos fotogramas del 2006 (Vuelo PNOA), mientras que la imagen de la derecha muestra el solapamiento entre dos fotogramas del año 1956 (Vuelo Americano Serie-B). 18 El hecho de que el solapamiento de los fotogramas usados sea de un 60% no solo implica que éste esté algo por debajo de lo recomendado y por tanto es posible que se reduzca la precisión de los resultados a la hora de generar un modelo, sino que también dificulta la búsqueda de puntos de control, ya que como se ha mencionado anteriormente, estos deben ser localizados en una franja donde coincidan al menos 3 fotogramas (Figura 8) y dicha franja será más estrecha cuanto menos solapamiento tengamos. Esto unido con el hecho de que la zona de estudio está conformada mayormente por un terreno con pocas edificaciones u otros elementos característicos con rasgos fácilmente identificables supuso un incremento en el grado de dificultad a la hora de encontrar puntos de control fiables. 3. Procedemos a buscar los puntos de control, es decir, puntos físicos localizados tanto en los fotogramas como en la ortofoto. Dichos puntos deben ser estables, es decir, que entre ambas fechas sus coordenadas XYZ, presumiblemente, no hayan variado a lo largo del tiempo (Figura 9). Estos puntos, pueden ser cualquier elemento natural o artificial sobre la superficie terrestre, como, por ejemplo, la esquina de un edificio, la base de un árbol, un cruce de caminos, etc., y tienen que estar claramente definidos y reconocibles en todas las fotografías donde aparezcan (Sanjosé et al., 2004). Para cada franja procuramos encontrar unos 2 o 3 puntos siempre y cuando sea posible; una vez localizados se digitalizan sobre la ortofoto de 2016, ya que es el mapa que estamos usando como referencia. De este modo, se acaba obteniendo una capa-inventario en formato vectorial con la posición de cada punto de control. En concreto se han identificado 18 puntos de control tanto para el año 2006 como para el 1956 y 19 para los fotogramas correspondientes a 2016. Figura 8. En esta figura se muestra en rojo el área de solapamiento entre dos fotogramas mientras que en amarillo se representa la franja de solapamiento entre 3 fotogramas, en la cual deben ser buscados los puntos de control. Figura 9. Ejemplo de un punto de control establecido, usando como punto físico estable la esquina de una alberca presente tanto en los fotogramas de 1956 así como en la ortofoto de 2016 usada como referencia. 19 4. El siguiente paso es asignar las respectivas coordenadas XYZ a los puntos de control. Esta tarea se realizará en un Sistema de Información Geográfica (ArcGIS). En el caso de las coordenadas X e Y el procedimiento es sencillo: mediante la herramienta “Add XY Coordinates” se selecciona la capa a la cual queremos asignar dichas coordenadas, y automáticamente al ejecutar la herramienta se crea en la tabla de atributos de la capa de puntos de control una columna con el valor de las coordenadas X y otra con el de las coordenadas Y de cada punto. 5. En cuanto a la coordenada Z, esta no es posible obtenerla a partir de la ortofoto ya que es un archivo en dos dimensiones, y por tanto solo contiene las coordenadas X e Y. Por esta razón, es necesario disponer de un archivo que contenga también valores de altitud. Es por ello que se utilizará la nube de puntos LiDAR 1ª Cobertura (2008-2015), ya que cada punto que conforma la nube sí contiene cada una de las 3 coordenadas. A partir de dichos puntos LiDAR, obtendremos un modelo digital de elevaciones (DEM), en otras palabras, un fichero ráster conformado por pixeles o celdas con valores de altitud. 4.2.1.1. Obtención de la coordenada Z. Antes de proceder con la explicación sobre cómo obtener la coordenada Z es importante tener en cuenta los siguientes aspectos sobre los archivos .las y .laz. El formato con extensión *.las es un formato estándar para el intercambio y distribución de nubes de puntos LiDAR. Posteriormente surgió el formato *.laz, que no es más que un formato comprimido del anterior, con el objetivo de reducir el tamaño de estos archivos. Precisamente el gran peso de estos archivos, debido a la alta densidad de puntos (habitualmente un mínimo de 0.5 puntos por m2), aconseja que la extensión geográfica que cubre cada uno sea bastante reducida, variando según el organismo que la distribuye. En nuestro estudio, los datos proceden del Instituto Geográfico Nacional (IGN), quien los distribuye en bloques de 2x2 km. Debido a esta limitación en la extensión de los archivos LAS/LAZ, suele ser habitual que para un área de estudio coincida con la intersección de dos o más bloques LiDAR (Orduña, 2019). Debido, por tanto, a la particularidad de los formatos de las nubes de puntos, es necesario para trabajar con este tipo de archivos disponer de herramientas adicionales, ya que por defecto en ArcGIS no están incorporadas. Por esta razón, se utilizó el paquete “LAStools”, que se descargó desde la web Rapidlasso GmbH (https://rapidlasso.com/), cuyas herramientas, en forma de archivos autoejecutables con extensión .exe, nos han permitido manipular y tratar tanto archivos *.las como *.laz (Orduña, 2019). Una vez obtenidas dichas herramientas y los archivos *.laz, descargados de la web del CNIG-IGN, se procedió a la elaboración del Modelo Digital de Elevaciones, siguiendo los pasos explicados a continuación: https://www.unigis.es/operaciones-lidar-lastools/ https://rapidlasso.com/ 20 1. Unir todos los archivos *.laz en un solo archivo mediante la herramienta “lasmerge”. Como ya se ha comentado previamente, los archivos descargados vienen en bloques de 2x2 km, por lo que es necesario unir todos ellos en uno solo, creándose un mosaico (Figura 10). 2. Para poder trabajar con el mosaico en ArcMap antes debemos descomprimirlo, obteniendo así un archivo en formato *.las. Para dicha operación usaremos la herramienta “laszip”. 3. El siguiente paso conlleva la creación de un fichero LAS Dataset en ArcGIS, es decir, un fichero que almacena la referencia a uno o más archivos *.las en nuestro disco para poder trabajar directamente con éstos sin necesidad de importarlos previamente en el entorno de ArcGIS. De este modo, podemos trabajar con estos archivos masivos de una forma más fácil y rápida (Tutorial ArcMap, 2020a). 4. Una vez disponemos del fichero LAS Dataset, se carga en dicho fichero el mosaico en formato *.las, creado anteriormente. Esta operación se realiza desde la pestaña “LAS Files”. Al incorporar nuestros archivos *.las al fichero LAS Dataset, este nos permitirá examinar dichos archivos en su formato nativo, ofreciéndonos estadísticas detalladas, clasificaciones de los puntos que conforman la nube y otros datos de interés. De estos, para el propósito de este estudio, nos interesa en especial el llamado “Point Spacing” o “Espaciado de los puntos”, que se necesitará posteriormente para la creación del MDE. Una vez hemos hecho esto, podemos cargar nuestro LAS Dataset en el entorno de ArcGIS, visualizar la nube de puntos y trabajar con ella. Gracias a las herramientas que nos proporciona ArcGIS podemos clasificar la nube de puntos según nuestro interés, así como aplicar filtros mediante los cuales podemos elegir que se muestren los puntos pertenecientes al primer retorno, los correspondientes al suelo, los no correspondientes al suelo o todos a la misma vez (Tutorial ArcMap, 2020a). Figura 10. Unión de los archivos *.laz en uno solo mediante la herramienta “lasmerge”. 21 5. Elaboración del Modelo Digital de Elevaciones (MDE). Mediante la caja de herramientas “LAS Dataset”, dentro del entorno de ArcMap, se emplea la funcionalidad de filtro y se selecciona la opción “Ground”, con el fin de utilizar solo aquellos puntos que se corresponden con el suelo desnudo o con el terreno (Figura 11). A continuación, con la herramienta “LAS Dataset to Raster” se procede a generar el MDE, estableciendo un tamaño de celda adecuado, en proporción a la densidad de puntos por metro cuadrado. En nuestro caso, ese tamaño sería de 4 m, ya que según se ha consultado en los manuales de ArcGIS se indica “que, aunque el espaciado de punto promedio puede parecer un tamaño de celda adecuado para el ráster de salida, este suele generar un archivo innecesariamente grande con demasiadas celdas vacías ya que los puntos LiDAR no se espacian de manera uniforme” (Tutorial ArcGIS Pro, 2020a). Por tanto, lo recomendable para garantizar un mejor resultado es utilizar un tamaño de celda que sea varias veces mayor que el espaciado de punto promedio, pero lo suficientemente pequeño como para identificar huecos resultantes al excluir aquellos puntos no correspondientes al suelo desnudo (árboles, edificaciones, etc.). Es por esto que se recomienda que el tamaño de celda sea de cuatro veces el espaciado de punto y, teniendo en cuento esto, en el área de estudio, el espaciado es aproximadamente de 1 m (0,953m). Por tanto, el tamaño de celda que debemos establecer sería de 4 m, esperando obtener un promedio de 16 puntos por celda. 6. Obtención de la coordenada Z. Una vez que se ha generado el MDE, el valor de la coordenada Z de la capa vectorial de puntos de control se extrae mediante la herramienta de ArcGIS “Extract Values to Points” perteneciente a la ToolBox “Spatial Analyst Tool” (Figura 12). Figura 11. Selección de los puntos pertenecientes al suelo desnudo mediante el filtro “Ground”. Figura 12. Imagen de la tabla de atributos de la capa vectorial correspondiente a los puntos de control una vez se han extraído las coordenadas X, Y y Z para cada uno de ellos. 22 4.3. Aerotriangulación o alineación de las fotografías aéreas o fotogramas. La siguiente etapa consiste en proceder a la alineación de las fotografías. Este proceso consiste en la estimación automática y precisa de la posición, rotación y propiedades de la cámara (distancia focal, punto principal y distorsión de la lente) para cada fotograma introducido (Bentley Systems, 2020a). Los puntos de control con sus coordenadas X,Y,Z (sistema de proyección ETRS89, UTM huso 30) se utilizarán, como ahora se explicará, para georreferenciar los modelos 3D de cada fecha. Los pasos para proceder a la alineación de las fotografías precisan la realización del siguiente conjunto de tareas: A) Importar los fotogramas desde el Menú de ContextCapture Master, y establecer si es posible los parámetros de la cámara fotográfica (“tamaño del sensor” y “distancia focal”). Siempre es recomendable especificar el valor de dichos parámetros, pero en caso de no disponer de ellos el software hará una estimación de estos de forma automática. Los parámetros de calibración de la cámara de los vuelos del PNOA están disponibles en el portal web del CNIG, aunque no especifican el tamaño del sensor de la cámara usada. También, para el vuelo de 1956/57 sólo se dispone de la distancia focal. Vuelo Modelo de cámara Distancia focal 1956/57 152.47 mm 2006 UltraCam D 105.20 mm 2016 UltraCam Eagle 100.50 mm Figura 13. Localización en ContextCapture de cada uno de los puntos de control en tres fotogramas distintos. Tabla 3. Parámetros de la cámara fotográfica, modelo de la misma y año del vuelo correspondiente. 23 B) Importar los puntos de control. Previamente, las coordenadas de los puntos de control de la tabla de atributos del archivo vectorial se han exportado a un fichero en formato de texto (con extensión .txt o .csv.). Una vez se tiene este fichero, se puede “cargar” directamente, desde ContextCapture, mediante la opción “Control points import wizard”. C) Localizar y establecer cada punto de control en sus correspondientes fotogramas. En la vista “Surveys” de ContextCapture, cada punto de control se posiciona en al menos tres fotogramas distintos (Figura 13). D) Una vez posicionados todos los puntos de control, el siguiente paso es ejecutar la aerotriangulación de las imágenes. En este paso es necesario arrancar el programa adicional “ContextCapture Engine”. Además, de forma manual se selecciona en las opciones de configuración una densidad de “keypoints” alta, es decir puntos de interés que detecta automáticamente el programa en las imágenes. Al terminar la aerotriangulación, se puede visualizar el modelo 3D formado por una nube de puntos comunes preliminar de baja densidad, que permite tener una primera aproximación del modelo 3D y la posición de las cámaras. En esta visualización se nos muestra además la posición y extensión que ocupa cada fotograma, así como la localización de cada Punto de Control (PC) (Figura 14). El color de los PC nos indica la calidad de la reproyección: verde (buena), amarilla (regular), roja (mala). Figura 14. Visualización de la posición de los fotogramas y los puntos de control sobre la nube de puntos 3D. 24 4.3.1. Control de calidad del proceso de aérotriangulación en ContextCapture. Para más información sobre la calidad de proceso de la aerotriangulación y del error medio cuadrático (RMSE, root mean square error) de los puntos de control se puede revisar el informe de calidad que se genera automáticamente (Bentley Systems, 2020b). Dicho informe nos muestra las principales propiedades y estadísticas de la aerotriangulación. Entre ellas podemos encontrar información sobre los resultados de la calibración de la cámara, representaciones visuales de la distorsión estimada de la lente o los errores de reproyección, etc. De toda la información que nos ofrece es de especial interés aquella relacionada con los puntos de control (PC) utilizados, ya que el reporte nos indica para cada PC el error medio cuadrático de la reproyección en pixeles, y los errores 3D, horizontal y vertical en metros. Mediante el análisis de estos datos podemos conocer cuál es la calidad y precisión de los puntos de control utilizados, y así poder detectar, en caso de que los hubiera, aquellos puntos con un error demasiado alto y que por tanto necesitan ser eliminados o reemplazados. En el caso del error medio cuadrático de la reproyección de los puntos, el programa nos indica de forma visual mediante colores el nivel de precisión de cada uno de ellos (Figura 15). Además, también será necesario conocer de cada modelo el error medio cuadrático global para el error 3D, ya que deberemos hacer uso de él más adelante a la hora de generar superficies de error, las cuales serán de gran importancia para la producción de los modelos digitales de diferencias de elevación (DoDs). Dicho error 3D es calculado a partir del error vertical y el error horizontal mediante la siguiente ecuación: (Error Vertical)2+(Error Horizontal)2= (Error 3D)2. Vuelo Error Horizontal Error vertical Error 3D 1956-57 X: 2,686; Y: 1,925 2,237 3,991 2006 X: 0,225; Y: 0,355 0.972 1.059 2016 X: 0,257; Y: 0,530 1,042 1,197 En nuestro caso de estudio se procuró no obtener ningún punto de control, cuyo error de reproyección excediera a más de 3 pixeles (punto en color rojo) (Figura 12) en ninguno de los modelos, aunque esto fue inevitable para los fotogramas de 1956-57. Se trata de Figura 15. Representación mediante colores del error medio cuadrático de la reproyección. Tabla 4. Errores medios cuadráticos globales de los errores 3D, horizontal y vertical de cada modelo. 25 fotogramas antiguos, de peor calidad. En este sentido, se presentan dos factores de importancia de los cuales depende la calidad de los productos obtenidos: el solapamiento de los fotogramas y la calidad de estos. Como ya se ha mencionado previamente, el solapamiento recomendado para una foto-reconstrucción es del 70%, mientras que el de nuestros fotogramas es aproximadamente del 60%, hecho que también limita la posibilidad de encontrar un mismo punto en tres fotogramas distintos o más y por tanto la fiabilidad de los puntos. En lo referente a la calidad de los fotogramas, determinada por su textura, que es la característica fundamental de la que depende la facilidad con la que es posible identificar elementos que puedan ser utilizados como puntos comunes. La calidad de imagen de los fotogramas históricos se ve afectada en muchos casos durante el proceso de conservación y de escaneado, el cual no siempre se realiza de una forma adecuada propiciando una pérdida de calidad en forma de zonas borrosas, manchas, presencia de artefactos, etc., y un incremento de errores en el proceso de alineación. En resumen, debido a estos factores la obtención de errores bajos de reproyección en los puntos de control obtenidos para fotogramas históricos resulta difícil (Llena et al., 2018). De esta forma, el error medio cuadrático (RMSE) global 3D fue, para los fotogramas de 1956, de casi 4 m (Tabla 4). En cambio, para los años 2006 y 2016, se obtuvo un RMSE global 3D más bajo, en torno a 1 m (Tabla 4). En los fotogramas de estos vuelos no se obtuvo ningún punto de control con un error de reproyección superior a 3 pixeles, (Figura 15), si bien, también es cierto que algunos de los puntos obtenidos fueron clasificados con un error de reproyección medio “amarillos” (Figura 15). 4.4. Producción del modelo 3D y de otros productos cartográficos. Una vez que se ha verificado que el RMSE de los puntos de control es aceptable en los tres vuelos, podemos proceder a presentar la producción de la malla 3D (3D mesh) o modelo 3D de cada vuelo. Previamente se establecieron los ajustes de la reconstrucción: por un lado, se optó por una precisión geométrica “Extra”, y, por otro, se eligió la opción “Adaptative Tilling” para la creación de las mallas de 2006 y 2016, debido a que el número de fotogramas en estos vuelos es mayor y ocupan bastantes megas (Tabla 5). Esto supone un aumento considerablemente del consumo de recursos computacionales, y puede, incluso, ocurrir que el ordenador no llegue a ser capaz de procesar el modelo. Por tanto, con la opción “Adapative tilling” el programa, en lugar de generar el modelo en un solo archivo, lo divide en varias secciones o “tiles”, procesando varios archivos. Cabe destacar que mientras los fotogramas del vuelo 1956-57 están en formato *.jpg, los fotogramas de 2006 y 2016 se presentan en formato *.ecw, característicos por su alta compresión. Vuelo Nº de fotogramas usados Peso total (MB) y tipo de formato de los fotogramas 1956-57 10 989 (*.jpg) 2006 16 427 (*.ecw) 2016 17 825(*.ecw) Tabla 5. Número de fotogramas usados en cada vuelo y peso total de los mismos. 26 Durante el proceso de producción de la malla o modelo 3D el programa analiza todas las imágenes en profundidad, con el fin de obtener el máximo número de puntos comunes posibles. Este proceso es el que más recursos consume, pudiendo tardar el procesado varias horas dependiendo del número de fotografías del vuelo. En este caso, las mallas generadas tardaron una media de 2 horas en producirse. Los resultados se pueden visualizar en 3D, directamente en ContextCapture, o hacer uso del software “ContextCapture Viewer”, que permite visualizar los productos 3D a mayor tamaño (Figura 16), e incluso en modo estéreo (Figura 17), y observar la profundidad del relieve con gafas anaglifo. Figura 16. Visualización del modelo 3D en ContextCapture Viewer. Figura 17. Visualización del modelo 3D en modo estéreo en ContextCapture Viewer. 27 3.4.1. Productos secundarios del procesamiento fotogramétrico. Una vez hemos generado el modelo 3D de cada vuelo, el siguiente paso es obtener, desde el mismo software (ContextCapture Master), distintos productos o materiales cartográficos secundarios de forma automática. Para ello simplemente, teniendo seleccionada la pestaña de nuestra reconstrucción, “clicamos” en el botón “Submit new production”, tras lo cual se nos abrirá una ventana desde la cual podremos elegir el producto que queremos generar mediante el botón “Purpose” (Figura 18). En este TFM, se han generado para cada modelo su correspondiente nube de puntos 3D en formato *.las con el fin de realizar una comparación entre ellas, la ortofoto en formato GeoTiff (*tiff), que nos será de ayuda a la hora de digitalizar los elementos geomorfológicos del área de estudio para cada año, y los MDS en formato Esri Ascii Grid (*asc), con los cuales podremos más tarde generar los modelos de diferencias de elevación o altura (DoD) con el fin de analizar los cambios producidos sobre el terreno entre los MDS de dos años distintos. Hay que tener en cuenta que tanto para la ortofoto como para MDS el programa no genera una única imagen, sino que el resultado se divide en un mosaico de imágenes, las cuales necesitaremos unir en una sola. Para ello se utiliza la herramienta “Mosaic To New Raster” que nos proporciona ArcGIS. Las resoluciones de los distintos productos, así como otras características pueden verse a continuación en las siguientes tablas. Destaca la gran cantidad de puntos que detecta el programa para la misma superficie en los tres vuelos (Tabla 6). Se observa, no obstante, que el número total de puntos que obtiene el programa de las fotografías aéreas históricas es bastante menor que los que se obtienen de las fotografías del PNOA de 2006 y 2016, cuyo número es aproximado. Como se dijo, las fotografías áreas históricas Figura 18. Procedimiento para la generación de productos secundarios en ContextCapture. 28 son en blanco y negro y proceden del escaneado de los contactos o negativos, mientras que las fotografías del PNIOA son digitales. En cambio, la resolución espacial o lado de celda, que ofrece el programa por defecto, tanto para las ortofotos como para los MDS es bastante alta (Tablas 7 y 8), obteniéndose resoluciones inferiores al metro. Material / AÑO Formato Resolución Tratamiento/Software Ortofoto 1956 GeoTiff (*tiff) 0,66 m Fotogrametría / Bentley ContextCapture Ortofoto 2006 0,50 m Ortofoto 2016 0,48 m Material / AÑO Formato Resolución Tratamiento/Software DSM 1956 Esri Ascii Grid (*asc) 0,66 m Fotogrametría / Bentley ContextCapture DSM 2006 0,50 m DSM 2016 0,48 m Fecha Número de puntos Point Spacing 1956 92.276.361 0,626 2006 165.267.071 0,468 2016 180.776.091 0,448 Tabla 6. Número de puntos y espaciado de puntos para la nube de puntos 3D de cada modelo. Tabla 7. Características de las ortofotos creadas a partir de cada modelo. Tabla 8. Características de las modelos digitales de superficies creadas a partir de cada modelo. 29 4.5. Reconocimiento y digitalización de los principales elementos geomorfológicos sobre las ortofotos. A partir de las ortofotos de los diferentes años se puede realizar un análisis comparativo simplemente con la observación de las mismas (Figura 19); si bien en este TFM, además, las ortofotos se han utilizado para reconocer y digitalizar las principales unidades geomorfológicas, ya que, tal como señalan Llena et al., (2018) “las ortofotos u ortomosaicos son imágenes rectificadas geométrica y geográficamente que permiten el estudio planimétrico de unidades geomorfológicas de interés, así como analizar su evolución (Llena et al., 2018). Desde este punto de vista, las ortofotos han permitido la elaboración de mapas geomorfológicos sencillos, mediante la digitalización de las unidades geomorfológicas de la llanura de inundación del río Guadalfeo. Así, se han cartografiado en cada una de las fechas estudiadas formas como: canales fluviales, lecho menor y lecho mayor. También se han identificado áreas de deslizamientos. A partir de estos mapas geomorfológicos se ha realizado en un SIG un diagnóstico y análisis de los cambios topográficos y morfológicos producidos a lo largo de los años en el tramo seleccionado del río. Figura 19. Comparación visual de un tramo del rio para los distintos años mediante el uso de las ortofotos. 30 4.6. Método para la elaboración de un Modelo Digital de diferencias de Superficie con el software Geomorphic Change Detection (GCD). Tal como se ha explicado en la introducción, el acrónimo DoD (DEM of Differencing) hace referencia, en geomorfología, a un modelo digital de diferencias de elevación o aturas (en el caso de comparar modelos digitales de superficies), las cuales pueden ser resultado de la resta de varios MDE o MDS (MDS en nuestro caso) de distinta fecha. Esta sencilla operación se puede implementar en cualquier SIG, sin embargo, el resultado podría contener muchas imprecisiones si el usuario no tiene en cuenta los errores de los MDE o MDS. Por esta razón, este objetivo fundamental del TFM se resuelve mediante la utilización de un software gratuito, diseñado y especializado en la detección de cambios geomorfológicos, mediante la resta de modelos. Además, nos permitirá conocer un software nuevo, y valorar sus funcionalidades. Este software, llamado “Geomorphic Change Detection” (GCD, http://gcd.riverscapes.xyz/), ha sido desarrollado por “Riverscapes Consortium’s” (https://riverscapes.xyz/), si bien en sus orígenes fue desarrollado por Joe Wheaton (Utah State University Department of Watershed Sciences) y James Brasington (The University of Waikato, New Zealand), principalmente para la detección de cambios morfológicos en ríos de gravas, aunque también puede ser usado para la detección de cambios en otros tipos ambientes. Teniendo en cuenta que cada MDE o MDS presenta cierta incertidumbre al ser representado (por ejemplo, RMSE), la habilidad del software para detectar cambios es altamente dependiente de dichas imprecisiones inherentes a cada MDE o MDS individual. Por tanto, el problema fundamental al hacer la resta entre dos modelos, radica en detectar si los cambios totales producidos son realmente cambios topográficos del terreno o si dichos cambios están dentro de los márgenes de error de los modelos. Es por ello que el software GCD proporciona un conjunto de herramientas para cuantificar las imprecisiones mencionadas de forma independiente en cada MDS, y, posteriormente las tiene en cuenta a la hora de generar el modelo con los cambios de elevación (DoD) (http://gcd.riverscapes.xyz/). Las imprecisiones o errores de los MDS que se han utilizado en este TFM están recogidos de los informes de calidad de cada una de las foto-reconstrucciones realizadas. En concreto, el dato que usaremos será el error 3D, el cual se da en metros, para el RMS global (Tabla 4 y 9). ERROR 3D (m) AÑO 3.991 1956 1.059 2006 1.197 2016 Tabla 9. Error 3D en metros para los MDS de cada año. http://gcd.riverscapes.xyz/ https://riverscapes.xyz/ http://gcd.riverscapes.xyz/ 31 El uso del software GCD puede hacerse desde una versión Standalone o, también, se puede implementar como una herramienta en ArcGIS. La principal diferencia reside en que mediante la implementación en ArcGIS podremos visualizar los resultados del análisis GCD tan pronto como sean generados, mientras que con la versión Standalone careceremos de un visor. Por tanto, el software, una vez descargado del portal web, se instaló en ArcGIS, como una funcionalidad más de la caja de herramientas. Desde ArcGIS, los pasos para realizar las restas de los DSM se explican a continuación, llevándose a cabo un total de 3 restas: • 2006 - 1956 • 2016 - 2006 • 2016 - 1956 1. Creación de un nuevo proyecto. En primer lugar, crearemos un nuevo proyecto desde la barra de herramientas de GCD en ArcGIS (Figura 20). 2. Preparación de los modelos digitales de superficies (MDS). Una vez creado el proyecto GCD, se puede importar los dos modelos que se vayan a restar. No obstante, antes de ello debemos tener en cuenta una serie de consideraciones y requerimientos bastante importantes en cuanto a formato, resolución y extensión de los MDS, ya que serán indispensables para que la resta se realice de forma correcta. De lo contrario los pixeles de los modelos se estarían restando de manera dispar. Estos requerimientos que tienen que cumplir las capas (en este caso, MDS) para poder operar con ellos en el programa GCD, son los siguientes: • El primer requerimiento que exige el programa es que los modelos introducidos estén en formato *.tif, ya que es el formato con el que trabaja la herramienta. Los MDS se generaron en formato Esri Ascii Grid (*asc), por lo que fue necesario exportarlos a formato *.tif. Figura 20. Creación de un nuevo proyecto GCD. 32 • El segundo requerimiento es que los modelos introducidos deben tener la misma resolución de píxel. En nuestro estudio, la resolución de los DSM se dejó por defecto (tabla 8). Por lo tanto, para solucionar este problema se realizó un reajuste de la resolución de forma sencilla mediante la herramienta “Resample” de ArcGIS, y, finalmente se estableció una resolución de píxel de un metro (Figura 21). • Finalmente, un tercer requerimiento, es que los modelos han de poseer exactamente la misma extensión espacial, debiendo sus límites estar conformados por números enteros. También es fundamental que dentro de esta extensión no entren zonas de los modelos con pixeles nulos o sin valor. Para ello se procedió de la siguiente forma: o Se digitalizó un nuevo shapefile poligonal con el área de estudio común a los tres DSM, y cuidando que no incluyera áreas con celdas sin datos. o Una vez dibujado el polígono, mediante la herramienta “Edit Sketch Properties”, se modificó las coordenadas X e Y de las esquinas del área de estudio, de modo que estuvieran formadas por números enteros (Figura 22). o Así, una vez tenemos el polígono, que define el área de interés, con las dimensiones deseadas se guardó, y se usó para recortar los MDS mediante la herramienta de ArcGIS “Clip” (Figura 23). Figura 21. Uso de la herramienta “Resample” en ArcGIS. Figura 22. Modificación de las coordenadas X e Y de la extensión para el shapefile del área de estudio mediante la herramienta “Edit Sketch Properties”. 33 3. Añadir los modelos digitales de superficie (MDS). Teniendo todos estos requisitos solventados, ahora sí, es posible proseguir con las herramientas del programa GCD. Al añadir los MDS en la tabla de contenidos de ArcGIS, las funcionalidades del programa GCD, ofrecen de manera automática una representación visual de cada MDS, que incluye el mapa hipsométrico coloreado, el mapa de sombreado (hillshade), y el mapa de la hipsometría y sombreado superpuestos, mostrando una representación bastante clara del modelo digital de superficies (Figura 24). Figura 23. Uso de la herramienta “Clip” en ArcGIS para recortarlos MDS a partir del shapefile creado. Figura 24. Distintas representaciones visuales del MDS ofrecidas por GCD al importar nuestro MDS. 34 4. Creación de la superficie de error para cada MDS. La superficie de error es un fichero ráster que representan la incertidumbre en el modelo y se aplica de manera uniforme sobre toda su extensión (GCD, en-línea). En este TFM, se ha considerado, como incertidumbre de los DSM, el error RMSE global 3D de los puntos de control utilizados para generar cada DSM, y cuyo valor es proporcionado por el informe de calidad de ContextCapture. Se puede considerar este error como el nivel mínimo de detección de cambios (Llena et al. 2018). Por tanto, a partir del valor RMSE global 3D de cada modelo, se generan sus correspondientes superficies de error. Así, cuando se procede a realizar la resta entre dos MDS, el programa GCD pide al usuario que seleccione los dos MDS y la superficie de error de cada uno de ellos. La superficie de error individual, teniendo en cuenta, por tanto, el error RMSE global 3D de cada modelo, se crea en la ventana del programa GCD “Create Error Surface For Entire DEM Extent”, donde se incluye en el campo “Uniform error value (m)”, el valor del error 3D obtenido en cada modelo (Figura 25). 5. Como último paso antes de proceder con la generación del DoD añadiremos una máscara del área de interés, de modo que el resultado será generado para dicha área en concreto. Para ello, uniremos en ArcGIS mediante la herramienta “Merge” los polígonos digitalizados de las llanuras de inundación de los dos años que queramos estudiar. Una vez hecho esto, añadimos el polígono resultante de la unión como una máscara dentro de la interfaz de GCD (Figura 26). Figura 25. Superficie de error creada para el MDS de 2006 usando el error 3D. Figura 26. Adición de máscara en GCD. 35 6. Generación de los modelos digitales de diferencias de alturas. Finalmente, el último paso es la creación propiamente dicha del “DoD” o Superficie de Diferencias. A la hora de elegir un método de análisis de las incertidumbres se ha establecido un umbral probabilístico de 0.95, de forma que el 95% de los datos (los cuales serán estadísticamente significativos) serán incluidos en la evaluación mientras que el resto será descartado. El resultado es el Modelo Digital de Superficie de Diferencias de Elevación o Alturas, que se muestra en la figura 27. 7. Información estadística que ofrece el programa GCD. Además del mapa, como resultado, GCD nos ofrece información estadística sobre los cambios producidos en términos de área, volumen y promedios verticales; que se expresan tanto de forma numérica, como gráficamente en forma de histogramas (Figura 28). Esta síntesis estadística sobre erosión/sedimentación es muy útil, y en un SIG no se muestra de manera tan directa, sino que habría que obtenerla indirectamente. Figura 27. Pasos finales para la creación del DoD y resultado obtenido (en este ejemplo entre 2006-1956. Figura 28. Información estadística de los resultados ofrecida por GCD. 36 En el siguiente capítulo se analizarán los resultados obtenidos de los modelos de diferencias de elevación, pero sólo se ciñeran a las unidades geomorfológicas (toda amplitud de la llanura de inundación y áreas de deslizamientos), que fueron digitalizadas a partir de las ortofotos. Por esta razón, los modelos de diferencias de elevación se recortaron, usando los archivos vectoriales generados para cada unidad geomorfológica, de modo que nos quedaremos únicamente con la superficie de los modelos donde se localizan dichas unidades. Más concretamente se han realizado las siguientes series de recortes, y posterior análisis de los resultados: - Llanura de inundación a partir de la resta entre los años 2006-1956 - Llanura de inundación a partir de la resta entre los años 2016-2006 - Deslizamientos a partir de la resta entre los años 2006-1956 - Deslizamientos a partir de la resta entre los años 2016-2006 37 5. RESULTADOS 5.1. Análisis geomorfológico de la evolución del río Guadalfeo a partir de la interpretación de las ortofotos. A partir del análisis visual en planta de las ortofotos se observa principalmente que el cauce varía notablemente de posición entre los distintos años. Mientras que en el año 1956 el cauce presenta una morfología algo más entrelazada y un mayor caudal de agua (aunque puede estar relacionado con la fecha concreta de cada fotografía), en los años 2006 y 2016 éste se muestra con una morfología meandriforme, propia de zonas donde la pendiente y la energía con la que fluye el río son bajas (Figura 29). Otro factor que podría tener influencia en la alteración de la dinámica fluvial en los años 2006 y 2016 es la construcción de presas en ciertos tramos del río (Figura 30), las cuales, al retener el flujo sedimentario pueden provocar cambios en el patrón fluvial, y favorecer que la lámina de agua fluya más bien por cauces únicos. Figura 29. Cartografía del cauce del río para los años 1956, 2006 y 2016. Figura 30. Recorte de una zona del tramo en la que se puede apreciar la existencia de una presa cuya construcción fue posterior a 1956. 38 Al tratarse a su vez de un río con un cierto carácter torrencial o de tipo rambla es también posible observar en todos los años canales múltiples o trenzados por los que no circula agua, especialmente en 2006 y 2016 en los tramos más abiertos del río (Figura 29 B). Dichos canales y su variación con los años son probablemente resultado de las avenidas producidas por lluvias torrenciales o por la fusión de la nieve de Sierra Nevada. En cuanto al lecho menor, constituido principalmente por depósitos que forman barras aluviales, es posible apreciar que con el paso del tiempo el cauce presenta importantes modificaciones laterales siendo está más notable en el año 2016. Así, se observa un desplazamiento lateral en determinadas zonas dando a su vez como resultado la formación de barras de gravas más estables o aparentemente no funcionales, las cuales han ido siendo progresivamente colonizadas por la vegetación, y propiciando un estrechamiento del lecho menor (Figura 31). Aunque también se observa, un desplazamiento del lecho menor en la ortofoto de 2016, sobre un sector de 1956 formado por el lecho mayor, y ocupado por la actividad humana (Figura 31 y 32). Por tanto, podemos observar tanto en la llanura de inundación áreas en las que dicho desplazamiento ha sido provocado por erosión del margen externo como zonas en las que este se ha producido por sedimentación en el margen interno (Figura 31) produciéndose, además, brazos muertos debido a estos procesos de desplazamiento. Por último, en referencia al lecho mayor, no se aprecia un gran cambio en cuanto al uso del suelo. Aunque es verdad que en algunas áreas con el tiempo han aparecido algunas construcciones más de las que ya existían en 1956, lo cierto es que dicho aumento no es excepcionalmente notable, por lo que la ocupación humana en el lecho mayor se ha mantenido más o menos igual a lo largo del tiempo. También es apreciable que en la mayor parte del lecho mayor se han mantenido las zonas de cultivos en aquellos tramos donde la erosión y la sedimentación lo han permitido (Figura 32). Figura 31. Desplazamiento lateral hacia la parte interna del lecho menor con establecimiento de barras colonizadas por vegetación. 39 Como resumen, en la Figura 33, vemos que, en relación a los cauces, la ortofoto de 1956 muestra el cauce más inestable, seguido por la de 2006 y siendo la ortofoto del año 2016, la que refleja un cauce más estabilizado y, aparentemente, más encajado en el lecho. Por otra parte, se ha obtenido, en el SIG, la superficie de las tres unidades geomorfológicas: lecho mayor, lecho menor y cauce o canal (Tabla 10). La evolución del lecho menor, formado por barras aluviales, muestra una reducción de su superficie en la ortofoto de 2016, con una mayor presencia de barras estabilizadas. El lecho mayor en las ortofotos de 1956 y 2006 muestra una superficie similar, ocupando aproximadamente el 63% del total de la llanura de inundación, mientras que en la ortofoto de 2016 su superficie aumenta, llegando a ocupar casi el 70% del total. De igual modo, en la ortofoto de 1956 la superficie ocupada por el cauce o lecho es mayor que en las ortofotos más recientes, lo que se refleja en un mayor caudal, lo cual parece lógico, teniendo en cuenta que la fecha del vuelo es del 18 de marzo de 1957. En cambio, en las ortofotos de 2006 y 2016 la superficie ocupada por el cauce es similar, y representa menos del 3% de la superficie total. Figura 32. Recorte en el cual se puede apreciar el mantenimiento de las zonas de cultivo a lo largo de los años excepto en aquellas zonas donde el lecho menor ha aumentado y ganado terreno al lecho mayor. 40 1956 2006 2016 Unidad Área ocupada (ha) % Área ocupada (ha) % Área ocupada (ha) % Lecho mayor 261,7 63,3 258,5 63.5 280 68,5 Lecho menor 133,4 32,3 138,7 34% 116,8 28,6 Cauce o Canal 18,6 4,5 9,9 2.4% 11,9 2.9 TOTAL 413,7 100 407,1 100 408,8 100 Figura 33. Cartografía de las principales unidades geomorfológicas fluviales para cada año. Tabla 10. Cálculo de la superficie ocupada por cada unidad geomorfológica para cada año. 41 5.2. Análisis de los cambios de elevación o alturas a partir de DoDs. A partir de los DoD producidos mediante el software GCD podemos identificar aquellas zonas en las que se ha producido erosión (valores negativos) y aquellas otras donde se ha producido una sedimentación (valores positivos). Sin embargo, también hay que tener en consideración que estamos trabajando con MDS, los cuales incluyen todos los elementos presentes en la superficie, tales como árboles, edificios, etc. 5.2.1 Análisis de los cambios de elevación o alturas en la llanura de inundación para los años 1956 y 2006. En primer lugar, el análisis del resultado estadístico, que el programa GCD ofrece en gráficos de tipo histograma, muestra, para este intervalo temporal, que los valores presentan una alta variabilidad y que están clasificados mediante 3 colores: rojo para los valores de diferencia negativos o dónde la altura se ha reducido, azul para los valores de diferencia positivos o donde la altura ha aumentado, y gris para aquellos valores que no son estadísticamente significativos y que por tanto han sido eliminados del análisis. En términos de superficie, en la figura 34 A, se observa que la mayoría de los valores se concentran a la izquierda, por tanto, se corresponden con valores negativos (o erosión), cuyo valor máximo está entorno al -21 m de cambio de elevación; de igual modo, el histograma que muestra el volumen de cambio, refleja este valor extremo como el de mayor pérdida de volumen (Figura 34 B). Por su parte, los valores positivos de superficie (o sedimentación) muestran dos picos o máximos, uno entorno a los 8 metros y otro entorno a los 24 metros de cambio de elevación, mientras que para el histograma de volumen solo presentan un pico de acumulación entorno al 8. Asimismo, los gráficos de barras correspondientes a los cambios de superficie y de volumen (Figura 34 A y B) también reflejan que la mayor parte de la superficie muestra Figura 34. Histogramas de área (A), volumen (B) con sus respectivos gráficos de barras y gráfico de barras de los valores medios de perdida y ganancia de altura (C), para el periodo 1956-2006. 42 un descenso o hundimiento, es decir una pérdida de altura. Esto se traduce, como se aprecia en la figura 34 C, en una profundidad media de perdida de terreno más elevada que la altura media de ganancia. La representación de los valores negativos y positivos que se muestran en los histogramas y diagramas de barras se expresan numéricamente en la Tabla 11. Así, entre 1956 y 2006 se aprecia, el predominio de los valores negativos, de modo que entre ambas fechas la pérdida de altura del terreno (o erosión) ha sido bastante mayor, lo que ha supuesto un volumen de terreno perdido superior a 12.000 m3. Cabe destacar que los valores de la Tabla 11 no son representativos del total del área de estudio, sino que en el cálculo solo se han tenido en cuenta aquellos pixeles con valores significativos dentro del umbral probabilístico que se estableció para crear la superficie de error, previa al cálculo del Modelo de Diferencias de Elevación (DoD 2006-1956), cuyos umbrales fueron de 3.991 m (RMSE en 1956) y 1.059 (RMSE en 2006) respectivamente. Estos valores se plasman cartográficamente en e el mapa de diferencias de elevación entre el DSM de 1956 y el DSM de 2006 (Figura 35). Los negativos, relacionados con procesos de erosión, se concentran en la margen inferior izquierda del tramo de río seleccionado, mientras que los valores positivos (o sedimentación) se encuentran más dispersos hacia la parte central y el margen derecho. El análisis en detalle del mapa de diferencias de elevación (DoD 2006-1956) y su comparación con las ortofotos de ambos años permite interpretar las posibles causas de los cambios geomorfológicos. Así, se aprecia que: a) Los valores positivos están relacionados prácticamente en su totalidad con zonas en las que ha crecido vegetación (Figura 36) y en algunas otras por acumulación de materiales por la acción antrópica (Figura 38). b) Los valores negativos están relacionados en parte con un tramo del rio en el que uno de los meandros parece migrar hacia el norte produciendo una clara erosión en el lecho mayor y consecuentemente un ensanchamiento del lecho menor. Si bien es cierto que esto es así, también hay algunas zonas con valores negativos pertenecientes al lecho mayor en las que no se observa ningún cambio evidente a simple vista donde persiste la presencia de zonas de cultivo (Figura 37). Área total con pérdida de altura (m2) 574.856 Volumen total de terreno perdido (m3) 12.116.521 Profundidad media de la perdida de terreno (m) 21,07 Área total con ganancia de altura (m2) 158.550 Volumen total de terreno ganado (m3) 2.049.157 Altura media de la ganancia de terreno (m) 12,92 Área total (m2) 733.406 Volumen total de diferencia (m3) 14.165.678 Tabla 11. Valores de área y volumen obtenidos para el DoD 2006-1956, así como de las medias de los valores positivos y negativos. Dichas sumas de valores sólo incluyen los datos dentro del umbral probabilístico. 43 Figura 37. Valores negativos del DoD correspondientes a una zona del cauce en la que se ha producido una aparente erosión. Figura 36. Valores positivos del DoD correspondientes a zonas en las que se ha producido un crecimiento de la vegetación Figura 35. Mapa de diferencias de elevación para la llanura aluvial entre los años 1956 y 2006 (Se muestra el MDS sombreado de 2006 como referencia). 44 5.2.2 Análisis de los cambios de elevación o alturas en la llanura de inundación para los años 2006 y 2016. El análisis, en primer lugar, de los histogramas y gráficos de barras obtenidos de la comparación del MDS-2016 y MDS-2006, vemos que los valores presentan una distribución mucho más centrada y un menor rango de variabilidad. Los valores, tanto negativos como positivos se concentran en un rango más bajo y cercanos al 0. Si, Figura 39. Histogramas de área (A), volumen (B) con sus respectivos gráficos de barras y gráfico de barras de los valores medios de perdida y ganancia de altura (C). Figura 38. Talud producido por acumulación de materiales por la acción antrópica. 45 observamos los histogramas y gráficos de barras correspondientes a las figuras 39 A y B, podemos ver que siguen predominando los valores negativos (o erosión) tanto en área como en volumen. Asimismo, la figura 39 C refleja que los valores medios de las diferencias negativas o positivas (profundidad o elevación respectivamente) son bastante más bajos y similares con respecto a los obtenidos para el DoD de 1956-2006, siendo estos de -4,72 metros para los valores negativos y 3,66 para los valores positivos. Los resultados observados mediante los histogramas y gráficos de barras se pueden observar en la Tabla 12 de forma numérica. Así, se pone de manifiesto que el área total con pérdida de altura representa el 85,4% con respecto al total, que se traduce en un volumen total de terreno perdido de cerca de 8.827.039 m3. La representación cartográfica del modelo de diferencias de elevación MDS-2016 y MDS-2006 (DoD 2016-2006) refleja claramente como los valores negativos, relacionados con procesos de erosión, predominan prácticamente sobre toda la extensión de la llanura de inundación, especialmente en la parte central y el margen derecho del mapa. En cuanto a los valores positivos, éstos se localizan, sobre todo, pequeños sectores del tramo central, tanto en la margen izquierda como derecha de la llanura, y también en el tramo final (Figura 40). Cabe destacar, que a diferencia del DoD 2006-1956 donde la mayor superficie de error afectaba a una amplia superficie de la llanura de inundación, en cambio en el DoD 2016-2006 la superficie de error es menor (RMSE de 1,197 y 1,059 m respectivamente), de modo es posible detectar con menos incertidumbre los cambios geomorfológicos que han ocurrido entre ambas fechas. De este modo, el análisis en detalle del mapa y comparándolo con las ortofotos se pueden mostrar los siguientes ejemplos de modificaciones topográficas de la llanura de inundación: a) El área representada con un color “azul intenso”, ubicada en la parte central de la llanura, corresponde al crecimiento de una línea de árboles (Figura 42). b) En cambio, la zona “azul claro” del tramo inferior de la llanura, en su extremo izquierdo, parece corresponder a la sedimentación provocada por el río sobre barras aluviales internas, y en ciertas zonas a un crecimiento por parte de la vegetación (Figura 41). Área total con pérdida de altura (m2) 1.870.895 Volumen total de terreno perdido (m3) 8.827.039 Profundidad media de las diferencias negativas (m) 4,72 Área total con ganancia de altura (m2) 319.870 Volumen total de terreno ganado (m3) 1.171.084 Altura media de las diferencias positivas (m) 3,66 Area total (m2) 2.190.765 Volumen total de diferencia (m3) 9.998.123 Tabla 12. Valores de área y volumen obtenidos para el DoD 2016-2006, así como de las medias de los valores positivos y negativos. Dichas sumas de valores sólo incluyen los datos dentro del umbral probabilístico. 46 c) Las zonas simbolizadas con “color rojo intenso” se corresponden con la perdida de espesor de vegetación o con la desaparición de esta (Figura 42) debido a procesos de erosión. Se observa la desaparición de vegetación a costa de la extensión del lecho menor. d) Las zonas donde los valores negativos son intermedios se corresponden en parte con zonas dónde se ha producido un encaje del cauce, así como con zonas donde se ha dado una cierta erosión en los lechos o una pérdida de vegetación (Figura 42). e) Las zonas de color “rojo claro” (o valores negativos menores) se corresponden con área de erosión generalizada pero menor, que se produce en el lecho menor y en el lecho mayor (Figura 42). Figura 40. Mapa de diferencias de elevación para la llanura aluvial entre los años 2006 y 2016 (Se muestra el MDS sombreado de 2006 como referencia). 47 5.2.3 Análisis de los cambios de elevación o alturas en los deslizamientos de la ladera del valle del río Guadalfeo para los años 1956 y 2006. De la misma manera que se ha procedido con el análisis de los cambios de elevación en la llanura de elevación, también se ha procedido a analizar los cambios producidos en una serie de deslizamientos que ocurren, principalmente, en la ladera de la margen derecha del valle del río Guadalfeo, para el periodo entre 2006 y 2016. Todos estos deslizamientos se encuentran en la parte oriental del área de estudio, y ocupan una superficie de 30,7 ha. El análisis de los histogramas y gráficos de barras generados, tras la resta del MDS- 1956 y MDS-2006, permite observar principalmente que, en todos ellos, gran parte de los valores son excluidos del análisis por no resultar estadísticamente signifcativos (color gris), y que se relaciona con el elevado RMSE que proporcionó la geooreferenciación de los fotogramas del vuelo americano. Figura 42. Los valores negativos más altos (rojo intenso) se corresponden con zonas con una pérdida de espesor de la vegetación mientras que el resto (rojo intermedio y claro) corresponden a zonas de erosión del lecho menor y mayor, así como a perdidas menores de espesor vegetativo. La zona de color azul (valores positivos) se identifica con el crecimiento de vegetación. Figura 41. Los valores negativos (rojos) representan la incisión provocada por el cauce tanto en el lecho menor, así como la erosión provocada en las barras aluviales mientras que los valores positivos (azules) se corresponden con zonas donde se ha producido sedimentación o un crecimiento de la vegetación 48 Tabla 13. Valores de área y volumen obtenidos para el DoD 2016-2006, así como de las medias de los valores positivos y negativos. Dichas sumas de valores sólo incluyen los datos dentro del umbral probabilístico. En todo caso, tanto para el área como para el volumen, se han detectado una mayor cantidad de cambios postivos (valores azules), los cuales se presentan con una alta variabilidad, mientras que los valores negativos (rojo) representan una proporción escasa en cuanto a su distribución espacial (Figura 43 A y B). En cuanto a la media de los valores verticales (negativos y positivos), vemos que estos son similares con valores cercanos a los 10 metros, aunque son algo superiores en el caso de los valores positivos (Figura C). En la tabla 13 se muestra la cuantificación tanto de los valores de superficie como los de volumen. En ambas variable, los cambios detectados son muy bajos. Debido a la superficie de error más alta del modelo de 1956 no podemos asegurar si la ausencia apenas de modificaciones geomorfologicas en este intervalo temporal (50 años) puede ser debido a una combinación, como se ha dicho, al alto error 3D (RMSE de 3.991 m) o bien a una cierta estabilidad del terreno. Área total con pérdida de altura (m2) 2.497 Volumen total de terreno perdido (m3) 22.266 Profundidad media de las diferencias negativas (m) 8,97 Área total con ganancia de altura (m2) 15.850 Volumen total de terreno ganado (m3) 190.939 Altura media de las diferencias positivas (m) 12,05 Area total (m2) 18.347 Volumen total de diferencia (m3) 213.205 Figura 43. Histogramas de área (A), volumen (B) con sus respectivos gráficos de barras y gráfico de barras de los valores medios de perdida y ganancia de altura (C). 49 La interpretación visual de la cartografia (Figura 44), sólo refleja modificaciones en varios deslizamiento durante el intervalo de tiempo entre los años 1956 y 2006. Así, en detalle, en el extremo norccidental se observan modificaciones positivas que apuntan a posible acumulación de material (Figura 45); mientras que en el borde suroriental se observa una pequeña área donde se detectan cambios negativos, en contacto con la llanura de inundación, y que podrian relacionarse con la zapa del río (Figura 46). Figura 45. Mapa de diferencias de elevación comparado con las ortofotos de 1956 y 2006. Figura 44. Mapa de diferencias de elevación para los deslizamientos entre los años 1956 y 2006 (Se muestra el MDS sombreado de 1956 como referencia). 50 5.2.4. Análisis de los cambios de elevación o alturas en los deslizamientos de la ladera del valle del río Guadalfeo para los años 2006 - 2016. En primer lugar, los histogramas y los gráficos de barras nos muestran para el conjunto de las áreas afectadas por deslizamientos, que los valores positivos predominan, reflejando una mayor variabilidad que los valores negativos (Figura 47). El histograma para el conjunto de la superficie afectada por deslizamientos presenta un pico de valores negativos en -1.5 metros mientras que para los valores positivos se observan dos picos; uno entorno a los 4 metros y otro entorno a los 9 metros (Figura 47 A). Para el histograma de volumen, los valores negativos siguen mostrando su máximo en el -1.5 así como para los valores positivos también se siguen observando ambos picos, aunque el pico entorno a los 4 aparece de forma más discreta (Figura 47 B). En cuanto a los valores medios de diferencia de elevación vemos que la media es mayor para los valores positivos mientras que para los valores negativos es mucho menor (Figura 47 C). Figura 46. Ampliación de la vista en un deslizamiento con valores negativos los cuales podrían relacionarse con la zapa del río (Se ha usado la ortofoto de 2006 como referencia). 51 La predominancia de los valores positivos es observable a su vez en la siguiente tabla que muestra los resultados estadísticos. La superficie de acumulación o ganancia es bastante superior en relación a la pérdida de terreno (Tabla 14). Por último, es posible observar de forma visual el predominio de las diferencias de elevación positivas en el DoD 2016-2016 generado. (Figura 48). La comparación de los resultados del DoD 2016-2006 con las ortofotos permite mostrar que en apariencia en las áreas afectadas por los deslizamientos no se ha producido un crecimiento de vegetación, aunque los resultados ofrecidos reconocen procesos de acumulación (valores positivos). Tan sólo se observa procesos de erosión en el deslizamiento más oriental (Figuras 49 y 50). Área total con pérdida de altura (m2) 26.925 Volumen total de terreno perdido (m3) 81216 Profundidad media de las diferencias negativas (m) 3,02 Área total con ganancia de altura (m2) 176.352 Volumen total de terreno ganado (m3) 1215747 Altura media de las diferencias positivas (m) 6,89 Area total (m2) 203.277 Volumen total de diferencia (m3) 1296963 Figura 47. Histogramas de área (A), volumen (B) con sus respectivos gráficos de barras y gráfico de barras de los valores medios de perdida y ganancia de altura (C). Tabla 14. Valores de área y volumen obtenidos para el DoD, así como de las medias de los valores positivos y negativos. Dichas sumas de valores sólo incluyen los datos dentro del umbral probabilístico. 52 Figura 49. Deslizamientos con valores azules positivos (sedimentación) predominantes. Figura 50. Deslizamiento con valores rojos negativos (erosión) predominantes. Figura 48. Mapa de diferencias de elevación para los deslizamientos entre los años 2006 y 2016 (Se muestra el MDS sombreado de 2006 como referencia). 53 6. CONCLUSIONES En este Trabajo Fin de Máster se ha mostrado la utilidad de la técnica de fotogrametría Structure from Motion (SfM) para producir modelos digitales de superficies y ortofotos de alta resolución a partir del procesamiento de fotografías aéreas históricas, cuyo tratamiento posterior nos ha permitido detectar cambios geomorfológicos en un tramo de la llanura de inundación del río Guadalfeo (Granada) y en áreas afectadas por deslizamientos. La aplicación de la técnica SfM a las fotografías aéreas del vuelo americano de 1956-57 nos ha permitido producir un modelo digital de superficies de una calidad bastante aceptable y apto para su uso a la hora de ser comparado con modelos digitales de superficies más recientes, si bien es cierto que el RMSE de la georreferenciación de dicho vuelo ha sido relativamente alto (cerca de 4 m), lo que ha dificultado la posibilidad de detectar cambios geomorfológicos fiables. En cambio, los MDS obtenidos a partir de los fotogramas digitales del PNOA (2006 y 2016) tienen un RMSE considerablemente menor, próximo a 1 metro. . Además, mediante la técnica SfM también ha permitido generar productos secundarios de gran calidad y utilidad, tales como las ortofotos de cada vuelo. Su uso y análisis ha permitido digitalizar las unidades geomorfológicas principales de la llanura de inundación (cauce, lecho menor y lecho mayor) y obtener y cuantificar su superficie. La interpretación visual de las ortofotos ha demostrado la enorme variabilidad del cauce del río Guadalfeo en los periodos analizados. La elaboración de Modelos Digitales de Diferencias de Elevaciones (DoD) se ha mostrado como un método eficaz para estudiar los cambios topográficos (superficie y volumen) de la llanura de inundación y de las áreas afectadas por los deslizamientos, ya demostrado en diversos estudios, como se indicó en la introducción. Para ello, las herramientas del programa GCD (Geomorphic Change Detection) simplifican el proceso de comparación de MDE o MDS y además expresan los resultados, tanto en gráficos (histogramas, gráficos de barras) como en tablas estadísticas y de manera cartográfica, que facilitan la interpretación de los procesos de erosión y sedimentación. El DoD 2006-1956 presenta una menor precisión que el generado para los años 2006 y 2016, pues ofrece valores de diferencia de alturas muy altos (tanto negativos como positivos), y excluye una gran parte del área de estudio, ya que la superficie de error es alta, ligado con el error de reproyección. En este sentido, los resultados tienen que ser analizados con cautela, Aun así, el mapa del DoD 2006-1956 permite detectar cambios de elevación, los cuales se pueden corroborar en parte mediante la comparación con las ortofotos, y se relacionan especialmente con el crecimiento o perdida de vegetación. En comparación, el DoD generado mediante la resta de los MDS de los años 2006 y 2016, nos proporciona una mayor cantidad de datos con una mejor precisión y una mejor identificación de los cambios producidos sobre el terreno, ya que los errores de reproyección eran menores en este caso. El mapa del DoD 2016-2006 refleja algunos cambios importantes relacionados con procesos de erosión. A modo de consideración final, es importante tener en consideración que el área de estudio en la que se ha realizado la foto-reconstrucción presenta una superficie grande, y estudiar la dinámica geomorfológica histórica precisa contar con fotografías aéreas, lo 54 cual supone trabajar con fotogramas de vuelos altos. Esto difiere si por ejemplo trabajáramos reconstruyendo un deslizamiento a escala local con fotogramas obtenidos mediante un vuelo con dron o con otras técnicas (TLS). Además, los cambios a estudiar en geomorfología fluvial, por lo general, son de una magnitud mucho menor que los estudiados en otros ámbitos como es el caso de los glaciares y, por tanto, el flujo de trabajo ha de seguirse rigurosamente para procurar minimizar todos aquellos errores relacionados con el proceso fotogramétrico (Llena et al., 2018). Así, la altura de vuelo y la calidad de los fotogramas históricos, junto con el bajo grado de solapamiento que estas poseen, dificulta la localización de puntos de control tanto por parte del usuario como por parte del software. Para una obtención de puntos de control de una manera más precisa que permitieran obtener un RMSE más bajo, lo más conveniente sería, realizar el apoyo de los fotogramas mediante mediciones con GPS diferencial sobre el terreno, estableciendo las coordenadas de forma precisa de aquellos elementos identificados en los fotogramas. Otra opción más sencilla, aunque menos eficaz, sería realizar una búsqueda de puntos de control mucho más exhaustiva sobre los fotogramas y la ortofoto de referencia. 55 7. BIBLIOGRAFÍA Asato, C.G., Perez Cerdán, F., Marín, G. (1996). SIG central del servicio geológico. La importancia del manejo de datos geológicos en formato digital. XIII Congreso geológico argentino. III Congreso de exploración de hidrocarburos, Capital Federal, Argentina. Ayala Carcedo, F.J., Corominas, J. (2003). Mapas de susceptibilidad a los movimientos de ladera con técnicas SIG. IGME. ISBN: 84-7840-466-X. 194 pp. Bentley Systems (2020a). ContextCapture User Guide. Dirección URL: https://docs.bentley.com/LiveContent/web/ContextCapture%20Help-v10/en/GUID- E71D0658-82CD-46F3-B400-17F95A3939EE.html (Fecha de consulta: 23-07-2020). Bentley Systems (2020b). ContextCapture Quality Report Manual. https://www.acute3d.com/QualityReportManual/en/v1.0/index.html (Fecha de consulta: 23-07-2020). Cook, K.L. (2017). An evaluation of the effectiveness of low-cost UAVs and structure from motion for geomorphic change detection. Geomorphology, , 278, 195-208. Coque, R. (1987). Geomorfología. Alianza Universidad textos.Madrid, 475 págs. Correa Muñoz, N.A. (2012). Método para la caracterización de las formas del terreno en zonas de montaña utilizando Modelos Digitales de Elevación. Caso: Departamento del Cauca. Tesis presentada como requisito parcial para optar al título de: Magister en Geomática. Díez Herrero, A.; Laín Huerta, L. y Llorente Isidro, M. (2006). Mapas de peligrosidad de avenidas e inundaciones, métodos, experiencias y aplicación. IGME. ISBN: 84- 7840-632-8. 229 pp. FAQS de Bentley. Bentley Communities. ContextCapture Frequently Asked Questions [en línea]. Dirección URL: (https://communities.bentley.com/products/3d_imaging_and_point_cloud_software/w/w iki/27678/contextcapture-frequently-asked-questions (Fecha de consulta:04-08-2020) Felicisimo, A.M. (1994). Modelos Digitales del Terreno. Introducción y Aplicaciones en las Ciencias Ambientales. Pentalfa, Oviedo, España. Fernández Chacón, F., Notti, D., Galve, J.P., Vicente Pérez; J. Azañón, J.M., Mateos, R.M., Lamas-Fernández, F., Monserrat, O., Roldán, F.J.; Pérez, J.L.; Colomo, C.M. y Gómez-López, J.M. (2015). Técnicas remotas para el análisis multiescala y multitemporal de fenómenos superficiales. XIV Reunión Nacional de Cuaternario, Granada. Fernández Lozano, J., y Guitiérrez Alonso, G. (2016). Aplicaciones geológicas de los drones. Revista de la Sociedad Gelógica de España 29 (1), 17. 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://docs.bentley.com/LiveContent/web/ContextCapture%20Help-v10/en/GUID-E71D0658-82CD-46F3-B400-17F95A3939EE.html https://docs.bentley.com/LiveContent/web/ContextCapture%20Help-v10/en/GUID-E71D0658-82CD-46F3-B400-17F95A3939EE.html https://www.acute3d.com/QualityReportManual/en/v1.0/index.html https://communities.bentley.com/products/3d_imaging_and_point_cloud_software/w/wiki/27678/contextcapture-frequently-asked-questions https://communities.bentley.com/products/3d_imaging_and_point_cloud_software/w/wiki/27678/contextcapture-frequently-asked-questions 56 Gomez, C. (2012). Historical 3D Topographic Reconstruction of the Iwaki Volcano using Structure from Motion from Uncalibrated Aerial Photographs. Hal-00765723, version 1 – 16 Dec 2012 GCD (en-línea). Geomorphic Change Detection Software. Riverscapes Consortium’s. http://gcd.riverscapes.xyz/Concepts/dod.html IGME (1979). Mapa Geológico de España. Escala 1:50:000. Hoja 1042 (Lanjarón). Instituto Geológico y Minero de España.Madrid, James, M.R., Robson, S. (2012). Straightforward reconstruction of 3D surfaces and topography with a camera: Accuracy and geoscience application. Journal of Geophysical Research: Earth Surface 117: F03017. doi: 10.1029/2011JF002289. Jiménez Perálvarez, J.D. (2005). Análisis de la susceptibilidad a los movimientos de ladera mediante un SIG en la cuenca vertiente al embalse de Rules, Granada. Memoria de Doctorado. Laín Huerta, L. (2002, ed.). Los sistemas de información geográfica en la gestión de los riesgos geológicos y el medio ambiente. Madrid, Instituto Geológico y Minero de España, pág. 288, ISBN 84-7840-458-9. Lane, S. N., Widdison, P. E., Thomas, R. E., Ashworth, P. J., Best, J. L., Lunt, I. A., Sambrook, G. H., Smith, Simpson, C. J. (2010). Quantification of braided river channel change using archival digital image analysis. Earth Surf. Process. Landforms 35, 971– 985. Llena, M., Vericat, D., y Martínez Casasnovas, J. (2018). Aplicación de algoritmos Structure from Motion (SfM) para el análisis histórico de cambios en la geomorfología fluvial. Cuaternario y geomorfología, 32 (1-2), 53-73. Micheletti, N., Chandler, J.H., Lane, S.N. (2015a). Structure from Motion (SfM) Photogrammetry. Geomorphological Techniques (British Society for Geomorphology), 2, 1-12. http://geomorphology.org.uk/sites/default/files/geom_tech_chapters/2.2.2_sfm.pdf Muñoz-Narciso, E., Béjar, M., Tena, A., Vericat, D., Ramos, E., Brasington, J., Gibbins, C.N. y Batalla, R.J. (2014). Generación de modelos topográficos a partir de fotogrametría digital automatizada en un rio de gravas altamente dinamico. XIII Reunión Nacional de Geomorfología, Cáceres, España. Nadal-Romero, E., Revuelto, J., Errea, P., López-Moreno, J.I. (2015). The application of terrestrial laser scanner and SfM photogrammetry in measuring erosion and deposition processes in two opposite slopes in a humid Badlands area (Central Spanish Pyrenees). SOIL 1561573 http://doi.org/10.5194/soil-1-561-2015. Nourbakhshbeidokhti, S., Kinoshita, A.M, Chin, A. and Florsheim, J.L. (2019). A Workflow to Estimate Topographic and Volumetric Changes and Errors in Channel Sedimentation after Disturbance. Remote Sensing 11 (5), 586. Orduña, F. (2019). Operaciones con datos LiDAR en LAStools. UNIGIS Girona, el programa de formación online en SIG del Servicio de Sistemas de Información http://geomorphology.org.uk/sites/default/files/geom_tech_chapters/2.2.2_sfm.pdf http://doi.org/10.5194/soil-1-561-2015 57 Geográfica y Teledetección (SIGTE) de la Universitat de Girona. Dirección URL: https://www.unigis.es/operaciones-lidar-lastools/ (fecha consulta: 26-07-2020). Rapidlasso GmbH (en linea). Creators of LAStools, LASzip, and PulseWaves. Blog at WordPress.com. Dirección URL: https://rapidlasso.com/ (Fecha de consulta: 23-07- 2020) Riverscapes Consortium’s (en-linea). Geomorphic Change Detection (GCD) Software. http://gcd.riverscapes.xyz/. Sanjosé Blasco, J.J., Martínez García, E., López González, M. (2004). Topografía para estudios de Grado. Bellisco, ediciones Técnicas y Científicas. Madrid, 413 págs. Smith, M.W., Carrivick, J. L., Quincey, D.J. (2016). Structure from motion photogrammetry in physical geography. Progress in Physical Geography, 40, 247-275. https://doi. org/10.1177/0309133315615805 Tomás, R. Riquelme, A., Cano, M., Abellán, A., y Jordá, L. (2016). Structure from Motion (SfM): una técnica fotogramétrica de bajo coste para la caracterización y monitoreo de macizos rocosos, 7. Tutorial ArcGIS Pro (2020a). Ayuda. Datos. Tipos de datos. Imágenes y ráster. Procesamiento y análisis. Funciones de conversión. Función De dataset LAS a ráster [en línea]. Dirección URL: https://pro.arcgis.com/es/pro-app/help/data/imagery/las-dataset- to-raster-function.htm (Fecha de consulta: 09-09-2020). Tutorial ArcMap (2020a). Administrar Datos. Tipos de datos. Dataset LAS. ¿Qué es un dataset LAS?. [en línea]. Dirección URL: https://desktop.arcgis.com/es/arcmap/10.3/manage-data/las-dataset/what-is-a-las- dataset-.htm (Fecha de consulta: 09-09-2020). Westoby, M.J., J. Brasington, J., Glasser, N.F., Hambrey, M.J. y Reynolds, J.M. (2012). ‘Structure-from-Motion’ photogrammetry: A low-cost, effective tool for geoscience applications. Geomorphology, 179, 300-314. Wheaton, J. M., J. Brasington, S. E. Darby, and D. Sear. (2010a). Accounting for uncertainty in DEMs from repeat topographic surveys: Improved sediment budgets, Earth Surf. Processes Landforms, 35, 136–156, doi:10.1002/esp.1886. Williams, R.D. (2012). DEMs of Difference. Geomorphological Techniques (British Society for Geomorphology), 2, 1-17. http://eprints. gla.ac.uk/114527/1/Williams%202012%20 DEMs%20of%20Difference.pdf. https://www.unigis.es/operaciones-lidar-lastools/ https://rapidlasso.com/ http://gcd.riverscapes.xyz/ https://pro.arcgis.com/es/pro-app/help/data/imagery/las-dataset-to-raster-function.htm https://pro.arcgis.com/es/pro-app/help/data/imagery/las-dataset-to-raster-function.htm https://desktop.arcgis.com/es/arcmap/10.3/manage-data/las-dataset/what-is-a-las-dataset-.htm https://desktop.arcgis.com/es/arcmap/10.3/manage-data/las-dataset/what-is-a-las-dataset-.htm