ANÁLISIS DE IMÁGENES PROCEDENTES DE VEHÍCULOS AÉREOS NO TRIPULADOS (UAV) PARA LA OBTENCIÓN DE APLICACIONES EN EL SEGUIMIENTO DE LA BIODIVERSIDAD. AUTORES: JOAN-CRISTIAN PADRO GARCIA ([email protected]) y JOAN TORRES BATLLÓ ([email protected]). DIRECTORES: Dr. LLUIS BROTONS ([email protected]) y MAGDA PLA ([email protected]). RESUMEN: El presente artículo trata la clasificación de las cubiertas del suelo en una zona quemada por el incendio del Solsonés (1998), a partir de imágenes obtenidas con UAV y un sensor con una resolución espectral centrada en el visible (RGB). A través de la realización de diferentes vuelos, el procesamiento de las imágenes capturadas y su correspondiente análisis, se ha realizado una clasificación de las cubiertas del suelo con una resolución de 25 cm de píxel corregistrada respecto a la ortofoto vigente a 25 cm de la cartografía oficial (ICGC) y una clasificación con diferenciación de especies vegetales de una resolución de 2.5 cm de píxel y un acierto del Índice Kappa del 0.9 como resultado final. El objetivo global del artículo consiste en proponer un protocolo solvente que marque el camino a seguir des de la adquisición de las imágenes hasta la generación de las clasificaciones de cubiertas del suelo, con tal de poder comparar las clasificaciones en años sucesivos y cuantificar la evolución de la regeneración post-incendio en la zona quemada por el incendio del 1998. PALABRAS CLAVE: UAV, fotogrametría, RGB, corrección geométrica, Modelos Digitales de Superficies, clasificación mixta. 2010), (Lelong et al, 2008). En general, la mayoría 1. INTRODUCCIÓN de autores que realizan clasificaciones mediante UAV han utilizado aparatos de tipo profesional La utilización de plataformas UAV en el marco como por ejemplo el ABS-Aerolight (Dunford et metodológico de la teledetección es al, 2009), Multi-rotor UAVs (Turner et al, 2014) i manifiestamente emergente, tal como apunta (Hung et al, 2014), MLB BAT3 UAV (Laliberte et Whitehead en su recuento de artículos académicos al, 2010). El hecho de utilizar UAV de un coste relacionados con esta tecnología publicados entre superior al nuestro les permite obtener imágenes el año 2000 y el 2013 (Whitehead et al, 2014). Las con menos error posicional (x-y-z). Los factores principales virtudes de los UAV respecto a las que marcan la diferencia son el sensor utilizado y el plataformas utilizadas en la teledetección, GPS de a bordo. En su caso, estamos hablando de típicamente aviones tripulados y satélites, caen en sensores como: Canon 550D Digital Single Lens la resolución espacial (menor a 2.5 cm) y la Reflex (Turner et al, 2009), Sony NEX 7 (Hung et resolución temporal (a la carta) que pueden al, 2014). Las características de la cámara utilizada ofrecer. En cambio, los principales inconvenientes en el artículo de Calvin Hung son (tabla 1): se encuentran en la autonomía del vuelo y la resolución espectral, ya que las plataformas clásicas Tabla 1: Características técnicas sensor Sony NEX tienen una capacidad de carga muy superior a las 7. de los UAV y por tanto pueden transportar Modelo Sony NEX 7 sensores mucho más complejos. En cuanto al Columnas y filas CCD 6000 x 4000 coste económico, las plataformas tradicionales Longitud focal (mm) 16 cobren mucha más superficie y por tanto el precio Campo de visión angular (º) 76 x 55 por hectárea levantada es más reducido; no Huella a 20 m (m) 30 x 2 obstante, los UAV tienen un coste de mercado Resolución a20m (mm/píxel) 5 asumible y eso los hace muy populares. Fuente: (Hung et al, 2014). Delante de las pocas experiencias que existen en la Para la adquisición y la planificación del vuelo cada generación de clasificaciones empleando autor, utiliza, en función de la plataforma escogida, únicamente el espectro visible, una de las un software que pertenece al fabricante de la aportaciones del artículo que tenéis en manos plataforma y con éste planifican el vuelo. Cada consiste en utilizar las regiones espectrales situadas plataforma trae incorporado un GPS de precisión entre los 0.4 micrómetros y los 0.7 micrómetros que ayudará a mosaicar las imágenes y a (RGB). La gran mayoría de estudios previos corregistrarlas. recurren a sensores capaces de registrar datos en el espectro infrarrojo (Feng et al, 2015), (Hung et al, Como el/la lector/a tendrá ocasión de observar en 2014), (Remondino et al, 2011), (Laliberte et al, el apartado material i métodos, las diferencias

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

cualitativas entre el tipo de plataformas y sensores utilizados por los autores de referencia y los utilizados en el artículo aquí introducido son considerables (tablas 3 i 4). La mayoría de autores utilizan Agisoft PhotoScan Pro ® 1.0.0 como el software apropiado para el procesamiento de las imágenes y la generación de ortofotos (Muñoz, 2014). Concluyen que calidadprecio es un programario muy asequible y da buenos resultados (Muñoz, 2014). Laliberte (Laliberte et al, 2008) utiliza como programa de fotogrametría digital eCognition. El factor determinante entre los autores de los artículos citados y este estudio está en las herramientas utilizadas (plataformas y sensores) y en el GPSdiferencial que llevan incorporados los UAV, hecho que les permite tener georreferenciadas las imágenes. También hay autores como Turner (Turner et al, 2008), Lelong (Lelong et al, 2008), Dandois (Dandois et al, 2013), Remondino (Remondino et al, 2011) y Zahawi (Zahawi et al, 2014) que utilizan GCP (puntos de control sobre el terreno) mediante GPS diferenciales, hecho que les permite obtener errores menores a la medida del píxel cuando mosaican las imágenes y corregistran las ortofotos. En este caso (Zahawi et al, 2014), utiliza un GPS-diferencial (Trimble GeoXT 2008).El error de este GPS es menor a 5 cm, hecho que permite obtener un MDE de muy buena calidad y una ortofoto georeferenciada. Otros autores (Remondino et al, 2011) y (Muñoz, 2014), exportan la nube de puntos y hacen una interpolación mediante ArcGIS.

Tabla 2: Variables utilizadas en la clasificación. RGB RATIOS (RGB) NDVI (RG,GB,BR) DESV_EST(R,G,B) VARIABLES (HSI) AREA DENSITY ROUGHNESS Fuente: (Laliberte et al, 2008). Una vez generadas estas variables el autor calcula la matriz de correlaciones para observar qué relación tienen entre ellas y discrimina todas aquellas superiores al 0.9 por medio del programa CART®. Acaba seleccionando un total de 10 variables de las cuales 6 son espectrales (NDVI, RG, Desviación Estándar G, Saturación y Ratio B, relación con R y Desviación Estándar B), 3 son espaciales (Área, Densidad y Curvatura) y 1 es característica contextual (diferencia de medianas de vecindad). Una vez realizada toda la clasificación evalúa la precisión mediante áreas de test y el Índice Kappa, extraído de la matriz de confusiones. El resultado del Índice Kappa para toda la zona de estudio fue de 0.61, con el que el autor se muestra satisfecho. El Centro Tecnológico y Forestal de Cataluña (CTFC) ha impulsado este proyecto en el cual se quiere evaluar qué posibilidades pueden ofrecer, en el ámbito de las aplicaciones ecológicas y forestales, la explotación de imágenes RGB procedentes de vehículos aéreos no tripulados (UAV). Una de las principales líneas de investigación del CTFC es el estudio de la regeneración post-incendio (Valor et al, 2015), (Tusell et al, 2014), (Castellnou et al, 2009), (Piqué et al, 2007), (Plana, 2007). El Trabajo pretende detectar las posibilidades de estas imágenes para el análisis de los procesos de regeneración forestal a partir de técnicas de fotogrametría y de la clasificación de cubiertas del suelo.

Las ortofotos generadas son verdaderas, es decir utilizan un Modelo Digital de Superficies. Teóricamente, este MDS puede ser comparado con MDS generados en posteriores vuelos para analizar las diferencias y extraer conclusiones sobre el crecimiento de la vegetación, así como analizar las diferencias con un MDE de detalle y estimar la cantidad de biomasa existente dentro del área de estudio. En cuanto a la clasificación en sí, Laliberte (Laliberte et al, 2008) realiza una clasificación jerárquica en función de las diferentes respuestas espectrales de los objetos/cubiertas que conforman su zona de estudio. El primer paso consiste en la separación de dos zonas: con sombras y sin sombras utilizando un Modelo Digital de Sombras. El segundo paso consiste en la diferenciación de diferentes sub-zonas dentro de la zona de no-sombra: vegetación y suelo desnudo.

Se define como objetivo general de este articulo el establecimiento de un protocolo que defina el camino a seguir des de la adquisición de las imágenes con UAV hasta la generación de clasificaciones de cubiertas del suelo y tipos de vegetación (figura 1) y que sea de utilidad en el CTFC. Un objetivo específico es contrastar cuales son los límites del espectro visible, hasta donde se puede llegar utilizando únicamente esta Fuente de datos para la realización de clasificaciones. Gracias al bajo coste del material utilizado, se pretende dar continuidad a los productos generados, realizando vuelos en años sucesivos para poder cuantificar (a través de un SIG) los cambios en los mapas de cubiertas del suelo generados. De este modo se espera conseguir información objetiva sobre la evolución de la regeneración post-incendio.

El tercer paso consiste en clasificar los diferentes tipos de vegetación, apoyándose en la creación de variables de textura que ayuden a discriminar cubiertas con respuesta espectral similar. El autor introduce las siguientes variables (tabla 2):

2

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

El área de estudio es una parte de la zona quemada en el incendio del Solsonés (1998). Éste, afectó un total de 26280 ha, de las cuales 15408 ha fueron bosques densos (Padró, 2014). La superficie estudiada ocupa 3 ha y se focaliza en un bosque situado en el municipio de Riner. Cabe remarcar que el área de estudio es meramente experimental, ya que el propósito de este artículo es mostrar un protocolo donde se definan los pasos a seguir en el tratamiento de imágenes obtenidas mediante UAV. Por tanto, no se analiza el ámbito de objeto de estudio en sí mismo, sino que se utiliza como muestra con la intención de poder repetir el mismo proceso en cualquier otro ámbito de características similares.

subsp. fagínea) con algún grupo de pinaza (Pínus nígra subsp. salzmánnii); la vegetación del sotobosque está formada brollas y matojos de media montaña. Una temperatura media anual de 13ºC y una pluviosidad anual 700 mm (Ninyerola et al, 2000) favorecen el establecimiento de estos bosques. Destaca la presencia de poblaciones de pinaza por encima de los rebrotes del roble. Geomorfológicamente, se encuentra justo sobre la divisoria de aguas de las concas del Cardener y el Segre, que en este tramo coincide con la carretera LV-3002. La altura del terreno se sitúa entre los 836 m y los 870 m (ICGC, full 374-642 LIDAR). Geológicamente forma parte de la Depresión del Ebro, concretamente la Formación Solsona, compuesta por materiales sedimentarios del oligoceno (Mata-Perelló, 2013).

Figura 1: Protocolo de trabajo general que se sigue en el presente artículo. Planificación de los vuelos

Figura 2: Localización y ámbito de la zona de estudio.

Adquisición de las imágenes

Procesamiento de las imágenes.

Generación de ortofotos

Corregistro entre vuelos y cartografía oficial

Generación y validación del MDS

Fuente: Ortofoto 25cm. ICGC y polígono de elaboración propia con MiraMon. Generación de variables y clasificación de cubiertas

3. MATERIAL Y MÈTODOS Fuente: Elaboración propia. El flujo de Trabajo general ha estado (figura 3):

En este trabajo, tanto la plataforma como el sensor son muy básicos (premisa del estudio), recurriendo a material de bajo coste para comprobar hasta qué punto puede explotarse y cuál es el máximo de provecho que se puede obtener.

a) Recopilación de la bibliografía con especial atención a los materiales y tipos de clasificaciones utilizados. b) Trabajo de campo basado en la planificación y realización de 4 vuelos, así como el inventario de especies vegetales presentes en la zona de estudio. La zona de estudio se ha visitado varias veces para observar los cambios biológicos que se producen debidos a la dinámica de los procesos naturales, así como las diferentes tonalidades que obtiene cada especie en función de la época y su comportamiento fenológico. La información recogida ha sido de mucha ayuda en el proceso de la clasificación (sección 3.3.4), aportando un criterio subjetivo que permite debatir los resultados (automáticos y objetivos) que el ordenador genera.

2. ÁREA DE ESTUDIO El área de estudio se localiza dentro de la comarca del Solsonés (Catalunya Central) situada entre la Depresión Central y el Prepirineo. Concretamente se trabaja en una zona del municipio de Riner, en el quilómetro 3 de la carretera LV-3002 (long: 1° 29' 55.48506'' lat: 41° 55' 22.92654'') y abasta un total de 3.3 ha (figura 2). La vegetación del estrato arbóreo que se encuentra en la zona de estudio, al igual que en gran parte del sud del Solsonés, está constituida mayoritariamente por bosques de roble valenciano (Quercus fagínea

c) Trabajo de estudio en el cual se procesan las imágenes obtenidas, se generan las ortofotos y los 3

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

MDT, se corregistran entre ellas y con la cartografía oficial, para finalmente realizar clasificaciones de las cubiertas del suelo.

-Sensor: El sensor utilizado es la cámara DJI FC 200. Sus características técnicas son (tabla 4): Tabla 4: Especificaciones técnicas del sensor. Modelo

Figura 3: Flujo del trabajo general. Recopilación bibliografía existente

Phantom Vision FC200

Trabajo de campo(generación del vuelo, observación de la zona de estudio, etc)

24

Fuente: Elaboración propia.

W: 110 N: 85

4384 x 3288

9 x 6.75

2.05292

Intervalos TemperaTiempo entre tura de de exposidisparos trabajo ción (s) (s) (ºC) 3 - 60

1/281

0 - 40

Fuente: DJI.

3.1 Material

- Imágenes utilizadas: Se han utilizado imágenes de 4 vuelos y la ortofoto oficial de 25 cm de resolución obtenida del ICGC. Las características generales de los vuelos son las siguientes (tabla 5):

Figura 4: Plataforma (DJI Phantom Vision 2 +) y sensor (DJI FC 200) utilizados.

Tabla 5: Características básicas de cada vuelo realizado. Vuelo Fecha Hora Nº imágenes totales local del vuelo 1 01/03/2015 14:24 84 2 03/03/2015 10:10 77 3 21/05/2015 10:00 134 4 25/06/2015 11:19 132 Fuente: Elaboración propia. Todas las imágenes obtenidas son en formato JPEG.

30 cm Fuente: Elaboración propia. -Plataforma: El UAV utilizado como plataforma para realizar los vuelos ha estado el DJI Phantom Vision 2 + (figura 4). Sus características técnicas son las siguientes (tabla 3):

Peso con batería (gr)

5

Profundidad de FOV (º) bits (b)

Trabajo de estudio(tratamiento de imágenes)

Tabla 3: Especificaciones plataforma utilizada.

Longitud Columnas DimenMedida focal y filas de siones de del píxel (mm) la CCD la CCD CCD (µm) (píxels) (mm)

técnicas

de

-Cartografía oficial: Hoja 273-100 de la Ortofoto de Catalunya 1:2500 (OF-25C) v3.3, vuelo realizado en Junio de 2014 y disponible en el sitio web del Instituto Cartográfico y Geológico de Cataluña (www.icgc.cat). Hoja 374-642 de la capa de puntos LIDAR disponible en el sitio web del Instituto Cartográfico y Geológico de Cataluña (www.icgc.cat).

la

Sensibilidad del planeamiento (m)

Máxima velocidad de giro horizontal (º/s) 1000 Vertical: 0.8 200 Horizontal: 2.5 Máxima Máxima velocidad Máxima longitud velocidad de de vuelo (m/s) (mm) ascenso (m/s) 6 10 350 Tiempo de Máximo peso de Temperatura de vuelo máximo despegue (Kg) trabajo (ºC) (min) 25 1.3 -5 a 50 Batería Peso de la batería Corriente de carga (gr) (A) 5200mAh Lipo 340 4 Frecuencia de Distancia de Sensibilidad del operación comunicación (m) receptor (Ghz) (1%PER) (dBm) 2.4 1000 -97 Alimentación Consumo de potencia 4 Baterías AA 100mA – 6V

-Software fotogrametría digital: El programario utilizado para realizar el tratamiento de las imágenes una vez obtenidas después del vuelo ha sido Agisoft PhotoScan. Las características del programa se encuentran en http://www.agisoft.com/ y en el manual del usuario (Agisoft, 2014). -Softwares GIS: Los programas utilizados para la realización del tratamiento de las ortofotos realizadas con Agisoft han sido MiraMon, ArcGIS, ENVI e IDRISI.

Fuente: DJI. 4

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

punto. Se estima un error en la determinación de las X, las Y y el conjunto de X-Y. Una vez hecho el proceso de la corrección geométrica y determinado el RMSxv, éste muestra un dato equivalente a la desviación típica. Este dato no ha de ser entendido como el error máximo que ha habido sino como un indicador del error global: asumiendo una distribución normal de los errores, el 68.27% de los errores se encuentran por debajo del RMS indicado, un 95.45% por debajo de 2 veces el RMS y un 99.76% por debajo de 3 veces el RMS. Generalmente y, para que este error sea aceptado, ha de estar por debajo de la medida del píxel (2 * RMS < medida píxel). En la evaluación de los Modelos Digitales de Superficies generados a partir de la interpolación de puntos, se ha utilizado la validación cruzada. Este método, también conocido como “jackknife” (Manly, 1986), realiza iterativamente la interpolación sin utilizar un punto (cambiante en cada iteración), midiendo la diferencia entre el valor conocido y el valor interpolado; la raíz cuadrada de la mediana de los errores al cuadrado, es el valor sintético de estimación del error en altitud (RMSz).

3.3. Metodología. Protocolo seguido Notas previas sobre las imágenes utilizadas y los resultados obtenidos: * El aspecto geométrico: Teniendo en cuenta el que el material disponible no incluye GNSS Diferenciales o GNSS RTK que determinen la posición exacta de ciertos puntos de control sobre el terreno con una precisión menor al lado de píxel más pequeño utilizado (2.5 cm), el mosaicado de las imágenes con las cuales se generará la ortofoto de cada vuelo se realiza de manera semi-automática con el programario de Agisoft. Éste, sitúa las imágenes individuales en el espacio basándose en la coordenada que el GPS de a bordo ha registrado en cada captura. No obstante, el GPS no tiene una precisión suficiente y una vez las imágenes mosaicadas y generada la ortofoto, ésta se encuentra desplazada respecto a la cartografía de referencia. La labor de corregistrar las imágenes entre ellas es complicada y no está libre de errores, ya que se fundamenta en el correcto posicionamiento de puntos de control digitalizados manualmente. Primero se han georreferenciado las ortofotos de alta resolución (2.5 cm) entre ellas tomando la del vuelo 3 como referencia; el interés cae en aprovechar al máximo la resolución que ofrecen los UAV. En segundo lugar, una vez densificadas a 25 cm, han sido todas georreferenciadas respecto la ortofoto oficial vigente a 25 cm, cartografía disponible en el sitio web del Instituto Cartográfico de Cataluña (www.ICGC.cat); el interés de la densificación cae en la obtención de ortofotos comparables con las oficiales en el espacio temporal que han entre vuelos del ICGC (típicamente 1 año), aprovechando el bajo coste y la versatilidad temporal de los UAV.

Respecto a la evaluación de los resultados vinculados con la clasificación temática se utiliza el Índice Kappa. Este índice evalúa la concordancia de medidas categóricas. En la interpretación del Índice Kappa (matriz de confusiones) se ha de tener en cuenta que este índice depende del acuerdo observado, pero también de la simetría del total estudiado. A través de la matriz de confusiones se recoge un acierto total de la clasificación y se ponderan los aciertos por el nombre de casos en cada categoría. La matriz de confusiones aporta información sobre errores de omisión y de comisión en cada categoría, y la exactitud del usuario y del productor. La exactitud del usuario cuantifica el porcentaje de celdas de test asignadas a una determinada categoría que realmente son de esta categoría, mientras que la exactitud del productor cuantifica el porcentaje de celdas de test iniciales de una determinada categoría y que han estado clasificadas correctamente como tales.

* El aspecto radiométrico: Cabe destacar que se trabaja con números digitales (DN), ya que se desconoce la función lineal que convierte los valores registrados por la CCD (sensor) de la cámara del UAV en radiancias; este hecho impide calcular la reflectividad de las cubiertas y las imágenes no son plenamente comparables; en definitiva, no se miden magnitudes físicas.

3.3.1. Adquisición de las imágenes La planificación del vuelo se ha realizado con el soporte del mismo programario que PHANTOM (DJI,2014) pone a disposición del usuario, en el cual se indica interactivamente el recorrido del vuelo deseado, la altura del vuelo deseado, la altura mediana respecto al terreno y el intervalo de disparos. Los parámetros escogidos han sido 5 m sobre el terreno y 3 segundos de intervalo entre captura y captura de imagen. También se ha probado un vuelo a 100 m pero ha sido descartado para mantener la serie a 50 m de todos los vuelos, que genera un tamaño de píxel teórico de unos 2

* Evaluación de la precisión: Para calcular la precisión en la georreferenciación de las imágenes, la interpolación de la nube de puntos y el resultado final de las clasificaciones se han empleado los siguientes métodos: Para determinar la precisión en la corrección geométrica de las imágenes, se han introducido puntos de test y se ha calculado la raíz cuadrada de la mediana de los errores (RMS) entre la posición conocida (a priori de la corrección) y la posición de destino (a posteriori de la corrección) de cada 5

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

cm (focal de 5 mm, escala 1:10000, medida de píxel en la CCD de 2.05 µm). Cabe destacar que el recorrido de la plataforma sigue en grandes rasgos las indicaciones establecidas a priori, pero no es una reproducción exacta, con la cual cosa la localización (x-y-z) del UAV en la captura de cada imagen se ha de obtener a partir del GPS de a bordo. Además, la muy variable actitud de la plataforma (cabeceo (φ), guiñada (κ) y balanceo (ω)) complica la geometría de captación de cada imagen. Como se ha visto en la introducción, estas variables se pueden saber con precisión utilizando tecnologías de precisión.

Figura 5: Localización de la captura de las 48 imágenes y estimaciones del error del vuelo 1:

3.3.2. Mosaicado de imágenes y generado de ortofotos. Agisoft genera un informe con el resumen de las localizaciones del UAV en la captura de cada imagen según la posición del GPS de a bordo. Este dato sirve para obtener unas localizaciones orientativas que el programa utiliza para ordenar las imágenes en el espacio. A partir de aquí, Agisoft orienta las fotos y busca puntos geométrica y radiométricamente parecidos que utiliza como puntos de enlace entre fotos individuales vecinas. A continuación se muestran las posiciones de captura de cada fotografía utilizada para hacer las ortofotos de los vuelos 1, 2, 3 y 4. Éstas, se han confeccionado con las fotografías centrales del área recubierta, ya que las más periféricas tenían muchas distorsiones, seleccionando por cada vuelo el número de imágenes pertinentes (tabla 6).

Fuente: AGISOFT. Figura 6: Localización de la captura de las 45 imágenes y estimaciones del error del vuelo 2:

Tabla 6: número de imágenes capturadas y orientadas para generar la ortofoto de cada vuelo. Vuelo Nº de imágenes Nº de imágenes totales del vuelo finalmente utilizadas 1 84 48 2 77 45 3 134 74 4 132 85 Fuente: Elaboración propia. Todas las imágenes obtenidas son en formato JPEG.

Fuente: AGISOFT. Figura 7: Localización de la captura de las 74 imágenes y estimaciones del error del vuelo 3:

En el informe que emite Agisoft figura la evaluación del error entre la localización que ha registrado el GPS de a bordo y la localización estimada que se utiliza para hacer el mosaico. En los siguientes gráficos, los errores en Z se representan con el color de la elipse, los errores en X-Y por la forma de la elipse (figuras 5 a 8).

Fuente: AGISOFT.

6

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

Figura 8: Localización de la captura de las 74 imágenes y estimaciones del error del vuelo 4:

Tabla 8: Características fotogramétricas vuelo 2. Númer Altura Resolució Área Puntos de o de de n sobre el levantad enlace imágene vuelo terreno a (ha.) establecid s median (m./píxel os por a sobre ) Agisoft el terreno (m.) 0.021087 45 51.359 2.81753 59362 1 Fuente: AGISOFT. Escala 1:10272 (focal 5 mm) Tabla 9: Características fotogramétricas vuelo 3. Númer Altura Resolució Área Puntos de o de de n sobre el levantad enlace imágene vuelo terreno a (ha.) establecid s median (m./píxel os por a sobre ) Agisoft el terreno (m.) 0.021433 74 52.203 5.63247 150069 9 Fuente: AGISOFT. Escala 1:10441 (focal 5 mm) Tabla 10: Características fotogramétricas vuelo 4. Númer Altura Resolució Área Puntos de o de de n sobre el levantad enlace imágene vuelo terreno a (ha.) establecid s median (m./píxel os por a sobre ) Agisoft el terreno (m.) 0.022099 85 53.824 6.08589 111234 4 Fuente: AGISOFT. Escala 1:10765 (focal 5 mm)

Fuente: AGISOFT. Los mosaicos de los dos primeros vuelos tienen importantes defectos de forma. La altitud del primer vuelo no ha sido la deseada (tabla 7), mientras que en el segundo vuelo el nombre de imágenes capturadas y útiles ha sido insuficiente para realizar un buen mosaicado (tabla 8). El vuelo 3 ha servido para realizar la mayoría de las experimentaciones debido a la suficiente cantidad de imágenes (tabla 9). No obstante, cuando se ha hecho el vuelo 4, se ha determinado como el más completo y útil para las finalidades de este trabajo (tabla 10). Este último vuelo no solo presenta un número de imágenes suficiente, sino que ha sido realizado en un día y hora en el cual la posición solar estaba muy cercana al zenit (tabla 5) y las sombras son prácticamente nulas, pudiendo calificar la imagen de limpia.

Máscaras: Una vez obtenidas todas las imágenes individuales seleccionadas, se aplica una máscara rectangular enmarcando la zona central de cada foto para eliminar las distorsiones de las zonas exteriores. Una máscara demasiado ancha reduce el solapamiento entre imágenes y entorpece la búsqueda automática de puntos de paso entre imágenes. Una máscara demasiado delgada no elimina las distorsiones de manera efectiva. La máscara utilizada ocupa 350 píxels en cada lado, es decir, se reduce prácticamente un 16% de los laterales (recubrimiento transversal) y un 21% de los márgenes superior e inferior (recubrimiento longitudinal); es por este motivo que el intervalo entre capturas es el mínimo posible, 3 segundos.

Tabla 7: Características fotogramétricas vuelo 1. Númer Altura Resolució Área Puntos de o de de n sobre el levantad enlace imágene vuelo terreno a (ha.) establecid s median (m./píxel os por a sobre ) Agisoft el terreno (m.) 0.017959 48 43.741 5.04906 85638 4 Fuente: AGISOFT. Escala 1:8748 (long. focal 5 mm)

Orientar fotos: Agisoft coloca las imágenes en el espacio siguiendo los datos registrados por el GPS del UAV en cada imagen en el momento de la captura, en el SRH lat/long y datum WGS84. Una vez situadas, automáticamente las orienta en función de los píxeles radiométrica y geométricamente homólogos que encuentra en las imágenes adyacentes. Este proceso es costoso en tiempo de procesamiento (>1hora). Una vez 7

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

orientadas, Agisoft proporciona un valor de ajuste entre fotos. En función de este valor y de un correcto recubrimiento del área a levantar, se seleccionan las imágenes más idóneas y vuelven a orientarse. Este paso es necesario para depurar las imágenes mal capturadas. Si hace falta, se repite el proceso varias veces. Exportar ortofoto: Se utiliza el formato TIFF para conservar la georreferenciación que ha realizado Agisfot. Aunque puede modificarse la dimensión del píxel de salida, se deja el píxel por defecto y des de MiraMon se controla el método de interpolación en el remuestreo. Se ha de ir con cuidado porque Agisoft no contempla que la CCD de DJI sea rectangular (como la utilizada en este trabajo), así que el lado definido en X tendrá que copiarse al de Y; en caso contrario el píxel de salida será rectangular. Importar ortofoto: Se visualiza la imagen con ENVI y se guarda en formato GeoTIFF con los datos intercalados por bandas (BSQ). Así, MiraMon es capaz de importarlo mejor al formato propio IMG, manteniendo la resolución radiométrica de 32 bits (8 x 3 bandas). Se utiliza siempre el valor 0 como a NODATA. Proyección cartográfica: Se reciben los ficheros con el SRH lat/long y datum WGS84, con el lado de celda por defecto (figuras 16 a 19). Estas serán las imágenes origen de los sucesivos procesos. Para facilitar los cálculos se opta por trabajar con un sistema de referencia horizontal UTM y datum ETRS89, así que se hace la proyección. Con el motivo de homogeneizar el lado de píxel, se trabajará con 2.5 cm (figura 9), ya que tiene la virtud de ser múltiple de los diferentes lados de píxel que queremos tratar posteriormente: 5 cm y 25 cm. Figura 9: Resultado de las ortofotos generadas por cada vuelo (2.5 cm de resolución).

Fuente: AGISOFT. Visualización con MiraMon.

8

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

Figura 12: Manca de corregistro entre dos ortofotos.

Por lo que refiere al método de interpolación para asignar el valor del píxel al hacer el cambio de proyección, hay 3 opciones: Vecino más cercano, bilineal (4 vecinos) o bicúbica (16 vecinos). En este caso se ha decidido utilizar el método de vecino más cercano, ya que la interpolación bicúbica genera valores irreales en los contornos de los objetos (figura 10), como por ejemplo una rama y el fondo sotobosque; en el mismo sentido, los valores NODATA en la imagen original que tienen algún vecino con valor válido, adquieren un nuevo valor, hecho no deseable. Como puede observarse, estos efectos no suceden con el método escogido (figura 11).

Fuente: Elaboración propia a partir de los vuelos 2 y 4. Corregistro (I): Se corregistran todos los vuelos respecto a la ortofoto del vuelo que presenta una mejor geometría. Se toma el vuelo 3 como referéncia y se busca el mismo píxel en las otras imágenes. Básicamente encontramos puntos en zonas geométricamente invariantes como las carreteras y sus líneas pintadas. Serán puntos de ajuste a los que se les asignará las coordenadas del píxel del vuelo 3. Se reserva un conjunto de puntos para evaluar la precisión del ajuste. Se utilizan el siguiente número de puntos de ajuste y de test, obteniendo los siguientes RMSxy (tabla 11):

Figura 10: Imagen original (izquierda) vs. Imagen georeferenciada con método bicúbico (derecha).

Tabla 11: Puntos de ajuste, de test i RMSxy de la corrección geométrica de la ortofoto de los vuelos 1, 2 y 4 respecto el vuelo 3. Vuelo Puntos Puntos RMSxy RMSxy ajuste de test (1er grado) (2on (m) grado) (m) 1 18 4 0.319440 0.1374 2 21 5 1.904840 0.507221 4 18 5 0.297586 0.099935 Fuente: Elaboración propia y MiraMon.

Fuente: Vuelo 3. Visualización con MiraMon. 2.5 cm/p Figura 11: Imagen original (izquierda) vs. Imagen georreferenciada con vecino más cercano (derecha).

La función utilizada para asignar a cada píxel orígen una celda destino puede ser de 1er o 2º grado. Se comparan los resultados de la georeferenciación mediante puntos de control, utilizando por un lado la función de 1er grado y por otro lado una función de 2º grado. La función de 1er grado únicamente hace inclinación, rotación y traslación mientras que la de 2º grado además de inclinar, mover y girar la imagen, también realiza cambios en la perspectiva. Es de esperar que la función de segundo grado ajuste mejor la imagen.

Fuente: Vuelo 3. Visualización con MiraMon. 2.5 cm/p Se ha de tener en cuenta, que aunque los RGB estén en el mismo SRH y lado de píxel, las imágenes no están corregistradas. La poca precisión del GPS de a bordo y la manca de puntos de control precisos sobre el terreno, hacen que haya un desplazamiento entre imágenes (figura 12).

FUNCIÓN DE 1ER GRADO: X’=A+Bx+Cy Y’=A’+B’x+C’y FUNCIÓN DE 2O GRADO: X’=A+Bx+Cy+Dx2+Ey2+Fxy Y’=A’+B’x+C’y+ D’x2+E’y2+F’xy Corrigiendo las imágenes con una ecuación de 1er grado se aprecia que en algunas zonas de la imagen el ajuste es excelente, pero en otras es peor. En cambio, con la ecuación de 2º grado el ajuste es visualmente mejor fijándonos en la alineación de 9

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

las carreteras. Además, cuantitativamente el cálculo del RMSxy da unos resultados mejores (tabla 11). Con tal de acotar el área de estudio, se recorta cada ortofoto con un marco, delimitando un ámbito general que excluya zonas periféricas. El marco con el cual se trabaja a partir de ahora tiene unas dimensiones de 175 x 175 metros (figura 13). Una vez corregistradas entre ellas, se pueden densificar las ortofotos de 2.5 cm a diferentes lados de píxel (5 cm, 10 cm, 25 cm), aunque el artículo solo presenta las de 2.5 cm y las de 25 cm. La série a 5 cm permitiría ofrecer un producto con imágenes más manejables sin sacrificar demasiada resolución espacial (tabla 12), mientras que la série a 25 cm tiene la finalidad de corregistrarla con la ortofoto del ICGC. Así, aunque perdiendo resolución espacial, las densificaciones de píxel tienen un interés concreto según los objetivos definidos.

ahora tomando como referéncia la ortofoto del ICGC a 25 cm de resolución. Se aconseja ajustar primero la imagen utilizada como referéncia (vuelo 3 en nuestro caso), y una vez ajustada aplicar los mismos coeficientes de la función de ajuste a las otras imágenes, en nuestro caso el vuelo 4. Con tal de escoger entre el uso de la función de 1er o 2º grado, y realizando las pruebas pertinentes se llega a la conclusión que en este caso se ajusta mejor la función de 1er grado (figura 14); basándonos en el criterio visual, se puede observar como las carreteras tienen mejor continuidad en la corrección con la función de 1er grado que con la de 2o grado (figura 15). Figura 14: Ajuste visual de la corrección geométrica del vuelo 3 respecto la ortofoto oficial utilizando una función de 1er grado.

Figura 13: Recorte donde se focaliza el ámbito de estudio.

Fuente: Vuelo 3 y ortofoto 25 cm ICGC. Figura 15: Ajuste visual de la corrección geométrica del vuelo 3 respecto la ortofoto oficial utilizando una función de 2o grado.

Fuente: Elaboración propia con la orto del vuelo 2. Tabla 12: Tamaño de cada fichero .IMG de las ortofotos (RGB) según el lado de celda. Lado de celda (cm) 2.5 5 25

Filas y columnas 7000 3500 700

Tamaño (Kbytes) 47.852 11.957 519

Fuente: Ortofotos recortadas por el marco comentado (175 x 175 metros) y tipo de dato Byte. Corregistro (II): En cuanto a las imágenes de 25 cm, ha de recordarse que las imágenes generadas hasta ahora no estan corregistradas con la cartografía de referéncia del ICGC. Por lo tanto, se ha de volver a realizar el proceso de georeferenciación con puntos de control pero

Fuente: Vuelo 3 y ortofoto 25 cm ICGC. 10

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

De hecho es una conslusión esperable, ya que solo se trata de trasladar la ortofoto des de una posición errónea a una posición correcta, no hace falta deformarla porque la proyección a UTM ya está hecha. Una vez se georeferencia la ortofoto del vuelo 3 se puede aplicar la misma función al vuelo 4. En los dos casos se han empleado los mismos 16 puntos de control y el RMSxy ha estado de 28.99381 cm. La corrección de los vuelos 3 y 4 a 25 cm con la cartografía de referéncia presenta un buen ajuste, ya que las ortofotos de los vuelos 3 y 4 estaban corregistradas con un error de 9.9 cm (tabla 11) y al aumentar el lado de píxel a 25 cm el error es menor a una celda en el 95% de los casos (figura 16). Figura 16: Detalle del ajuste visual entre la ortofoto a 25 cm vigente y el vuelo 3.

Figura 18: Detalle de la corrección geométrica por interpolación bicúbica del vuelo 3 (25 cm).

Fuente:Zoom a la ortofoto del vuelo 3 corregistrada con la ortofoto oficial con interpolación bicúbica. 3.3.3. Generación de Modelos Digitales del Terreno (MDT) Notas previas sobre los MDT, el objetivo de este apartado y los resultados obtenidos: Delante de la habitual confusión en el uso de la terminología empleada en este apartado, utilizamos los términos fundamentales tal como se definen en el “Diccionari terminològic de teledetecció” (Pons y Arcalís, 2012; entr. 1922, 1926 i 1927). En resumen:

Fuente: Elabrotación propia a partir de la ortofoto vigente a 25 cm ICGC y la ortofoto del vuelo 3 georeferenciada con puntos de control.

Modelo Digital del Terreno (MDT): Estructura

En cuanto al método de interpolación, visualmente llegamos a la conclusión que el método del vecino más cercano es el más apropiado, ya que mantiene más los contrastes (figura 17). En cambio, las interpolaciones lineal y bicúbica dan unos resultados más borrosos (figura 18). En el caso de un ráster con valores continuos como el del MDS, es preferible el método bicúbico para dar continuidad a la variable, pero tratándose de celdas que contienen valors RGB creemos que hace falta mantener la radiometria, aunque se trate de DNs (no de radiancias). Conviene tener presente que hasta este punto se han realizado dos procesos de interpolación (1 canvio de resolución + 1georeferenciación).

Modelo Digital de

numérica de datos que describe la distribución espacial de una variable cuantitativa continua.

Elevaciones (MDE):

Modelo Digital del Terreno que recoge los datos topográficos i/o batimétricos. Las alturas recogidas son ortométricas, no elipsoidales. Si no se dice lo contrario, se supone que los datos son sobre el suelo desnudo, sin vegetación ni construcciones humanas.

Modelo Digital de Superficies (MDS): Modelo Digital del Terreno en que la variable cuantitativa continua es la altura de la superficie, de manera que considera la vegetación, edificios, puentes, etc.

Agisoft genera una ortofoto verdadera, razón por la cual necesita un MDS. Este modelo de alturas lo genera automáticamente mediante trigonometría y paralaxi, con los mismos métodos que utiliza la estereoscopía (Agisoft, 2014), pero tiene una resolución entorno a los 30 cm y no se puede dominar el método de interpolación. Las alturas son relativas a la altura del vuelo, ni elipsoidales ni ortométricas. Es importante remarcar que aunque Agisoft denomina “DEM” a su modelo, se trata de un Modelo Digital de Superficies (MDS), es decir, no reproduce la altura del suelo sino de los elementos que contiene (árboles, edificaciones…). Siguiendo el artículo de Muñoz (Muñoz, 2014) se exporta la nube de puntos (ASCII (XYZ)) y se calcula un nuevo MDS con MiraMon.

Figura 17: Detalle de la corrección geométrica por vecino más cercano del vuelo 3 (25 cm).

Fuente: Zoom a la ortofoto del vuelo 3 corregistrada con la ortofoto oficial por vecino más cercano.

El objetivo de ese apartado consiste en evaluar la precisión del MDS generado a partir de la nube de 11

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

puntos emitido por Agisoft, y por lo tanto, ver si puede utilizarse con un error aceptable. La evaluación se lleva a cabo internamente mediante la validación cruzada como externamente mediante la comparación con un MDS de referencia. Todos los MDS presentados se han calculado a partir de la interpolación de la nube de puntos emitida por Agisoft, excepto el MDS de referencia que se ha calculado a partir de la nube de puntos LIDAR descargada del ICGC.

Figura 19: Sección de la hoja 374-642 LIDAR del ICGC y marco del área de estudio.

Un correcto ajuste del MDS generado a partir de la nube de puntos de Agisoft, tendría que permitir obtener los siguientes productos: -Estimación de la biomasa presente: restando un MDE de referencia con el MDS generado en un vuelo, teóricamente se puede calcular la altura de la vegetación sobre el terreno y por tanto derivar su volumen de biomasa.

Fuente: Hoja 374-642 LIDAR del ICGC. Figura 20: MDS a partir de los puntos LIDAR mediante interpolación por método IDW.

-Estimación del crecimiento de la vegetación: restando el MDS de un vuelo y el MDS de otro vuelo se podría estimar el aumento de biomasa y la localización espacial de la misma. -Cálculo del Modelo Digital de Iluminaciones: a partir de la posición solar en el momento de la captura de las imágenes, se puede asignar a cada píxel un valor que cuantifique el ángulo de incidencia de la luz solar respecto la superficie, información útil para poder corregir los efectos de la iluminación en una hipotética corrección radiométrica si finalmente se conoce la función lineal que convierte los DN en radiáncias. -Cálculo del Modelo Digital de Sombras: a partir de la posición solar en el momento de la captura de las imágenes, se puede deducir la sombra proyectada por los objetos de la superficie (en este caso la vegetación) y corregir así efectos radiométricos y topográficos no deseados como las zonas con sombras.

Fuente: Elaboración propia a partir de la hoja 374642 LIDAR del ICGC. Figura 21: MDS a partir de los puntos LIDAR mediante interpolación por método de esplines.

-Soporte a la clasificación: complementar con estos modelos las variables utilizadas para la clasificación de las cubiertas del suelo. Generación del MDS de referencia: Para obtener un Modelo Digital de Superficies de referencia, nos valemos de la base de puntos de LIDAR descargable en el sitio web del ICGC. Utilizaremos la hoja 374-642, la cual supera el ámbito deseado y por lo tanto recortaremos (figura 19); detectamos en este trabajo algunos problemas anómalos, ya que sabemos a priori que las alturas están alrededor de los 800 y 900 m.s.n.m. y hay algún punto con valor 275, 600…así que se depura el fichero con estos puntos anómalos (puntos (Z)DBF). • Interpolación manual de la nube de puntos. ("jackknife"). • Corregistro entre MDS para la comparativa entre fechas.

• Generación de variables (a partir del RGB, índices, ratios, texturas (5x5). • Generació de una matriz de correlaciones y de componentes principales para ver la correlación entre elles. • Generación de áreas test y áreas de entrenamiento para la clasificación. • Clasificación no supervisada con variables escogidas. Clasificación mixta con las áreas de entrenamiento. • Validación mediante áreas de test y toma de decisiones a partir de la matriz de confusiones.

Fuente: Elaboración propia.

23

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

BIBLIOGRAFÍA

monitoring of wheat crop in small plots”. Sensors, 8, 3557-3585.

AGISOFT. (2014). “Agisoft PhotoScan users manual: Professional Edition, 1.1”. http://downloads.agisoft.ru/pdf/photoscanpro_0_9_0_en.pdf

MANLY, B.F.J. (1986). Multivariate statistical methods. London: Chapman and Hall. MATA-PERELLO, J.M. (2013). “Recorregut de recerca geològica i minera per les comarques del Solsonès, La Segarra i La Noguera: Des de Solsona a Sant Climenç, Sanaüja i a Vilanova de l’Aguda”.

CASTELLNOU, M.; PIQUÉ, M.; MIRALLES, M.; PAGÊS, J. (2009). “Integració del risc d’incendi en la planificació forestal: règim d’incendis i incendis tipus a Catalunya”. Silvicultura, 60: 5-8.

MUÑOZ, E. (2014) “Aplicación de fotogrametria digital automatitzada en geomorfologia: dinàmica morfo-sedimentaria en ríos y gravas”. Projecte de final de Master en SIG i Teledetecció. Universidad de Zaragoza.

DANDOIS, J.; ELLIS, E. (2013). “High spatial resolution three-dimensional mapping of vegetation spectral Dynamics using computer vision”. Remote Sensing Environtment, 136, 259-276.

NINYEROLA, M.; PONS, X.; ROURE, J.M. (2000). "A methodological approach of climatological modelling of air temperature and precipitation through GIS techniques", International Journal of Climatology, 20, 1823-1841.

DJI (2015). “Phantom 2 Vision Plus User Manual v1.8 en V2”. http://download.djiinnovations.com/downloads/phantom-2vision/en/Phantom_2_Vision_User_Manual_v1.8 _en.pdf

PADRÓ, J.C. (2014). “Anàlisi dels canvis de les cobertes del sòl durant el període 1956-2009 en dues finques afectades per l’incendi de 1998 a la Catalunya Central”. Treball de fi de grau. UAB: Dipòsit Digital de Documents.

DUNFORD, R.; MICHEL, K.; GAGNAGE, M.; PIÉGAY, H.; TRÉMELO, M.L. (2015). “Potential and constraints of Unmanned Aerial Vehicle Technology for the characterization of Mediterranean riparian forest”. International Journal of Remote Sensing, 30:19, 4915-4935.

PIQUÉ, M.; ROMÀ, J.; VERICAT, P. (2007). “Caracterización de las masas de robles submediterráneos (Quercus humilis e híbridos afines) como base para la definición de propuestas de gestión en el prepirineo catalán”. Cuadernos de la Sociedad Española de Ciencias Forestales, 21, 67-74.

FENG, Q.; JIANTAO, L.; GONG, J. (2015). “UAV Remote sensing for urban vegetation mapping using random forest and texture analysis”. Remote Sensing, 7, 1074-1094.

PLANA, E. (Ed.) (2007). Gestió del risc d’incendi, ecologia del foc i restauració de zones cremades. Solsona: Xarxa ALINFO.

HUNG, C.; XU, Z.; SUKKARIEH, S. (2014). “Feature learning based aprroach for weed classification using high resolution aerial images from digital càmera mounted on a UAV”. Remote Sensing, 6, 12037-12054

PONS, X.; ARCALÍS, A. (2012). “Diccionari terminològic de Teledetecció”. Barcelona: Institut Cartogràfic de Catalunya.

HARALICK, R.; SHANMUGAM, K.; DINSTEIN, I. (1974). “Textural features for Imagen classification”. IEEE, 6 (3), 610-621.

PONS, X. (1996). "Estimación de la Radiación Solar a partir de modelos digitales de elevaciones. Propuesta metodológica". A: VII Coloquio de Geografía Cuantitativa, Sistemas de Información Geográfica y Teledetección. Juaristi, J. i Moro, I. (eds.) Vitoria-Gasteiz.

LALIBERTE, A.; HERRICK, J.; RANGO, A.; WINTERS, C. (2010). “Acwuisition, orthorectification, and Object-Based Classification of UAV imagery for rangeland monitoring”. Photogrammetric Engineering & Remote Sensing, 76 (6), 661-672.

REMONDINO, F.; BAZZARETTI, L.; NEX, F.; SCAIONI, M.; SARAZZI, D. (2011). “UAV photogrammetry for mapping and 3D modelling – current status and future perspectives”. International Archieves of the Photogrammetry, Remote Sensing and Spatial Informations Sciences. 38 (1)

LALIBERTE, A.; RANGO, A. (2011). “Imagen processin and classification procedures for analisys of sub-decimeter imagery acquired with an UAV over arid rangelands”. GIScience & Remote Sensing, 48 (1), 4-12.

TURNER, D.; LUCIEER, A.; MALENWSKÝ, Z.; KING, D.H.; ROBINSON, S.A. (2014). “Spatial co-registration of ultra-high resolution visible, multiespectral and thermal images acquired

LELONG, C.; BURGER, PH.; JUBELIN, G.; ROUX, B.; LABBÉ, S.; BARET, F. (2008). “Assesment of UAV imagery for quantitative 24

Análisis de imágenes procedentes de UAV para la obtención de aplicaciones en el seguimiento de la biodiversidad

with micro-UAV over Antartic Moss beds”. Remote Sensing, 6, 4003-4024. TUSELL, J.M.; VERICAT, P.; FRIGOLA, P. (2014). “Treballs de regeneració dels boscos cremats. El cas de la sureda”. A: TUSELL, J.M. i VERICAT, P. (eds.). XXXI Jornades Tècniques Silvícoles Emili Garolera. 64-73. VALOR, T.; GONZALEZ-OLABARRIA, J.R.; PIQUÉ, M. (2015). “Assessing the impact of prescribed burning on the growth of European pines”. Forest Ecology and Management, 343, 101-109. WHITEHEAD, K.; HUGENHOLTZ, C. (2014). “Remote sensing of the environament with small unmanned aircraft Systems (UASs), part 1: a review of progress and challenges”. J. Unmanned Veh. Systems, 7, 69-85. ZAHAWI, R.A.; DANDOIS, J.P.; , HOLL, K.D.; NADWODNY, D.; LEIGHTON, J.; ELLIS, C. (2014). “Using lightweight unmanned aerial vehicles to monitor tropical forest recovery”. Biological conservation, 186, 287-295.

CARTOGRAFIA OFICIAL UTILIZADA: INSTITUT CARTOGRÀFIC I GEOLÒGIC DE CATALUNYA. Full 273-100 de l’Ortofoto de Catalunya 1:2500 (OF-25C) v3.3, vuelo realitzat al juny de 2014 i disponible al lloc web de l’Institut Cartogràfic i Geològic de Catalunya (www.icgc.cat). INSTITUT CARTOGRÀFIC I GEOLÒGIC DE CATALUNYA. Full 374-642 de la capa de punts LIDAR disponible al lloc web de l’Institut Cartogràfic i Geològic de Catalunya (www.icgc.cat).

25