Aplicación de datos LiDAR a la gestión forestal

Grado en Ingeniería Geomática y Topográfica UNIVERSITAT POLITÈCNICA DE VALÈNCIA. Trabajo Fin de Grado en Ingeniería Geomática y Topográfica Alumno: Rafael Llorens Company Tutor: Prof. D. Alfonso Fernández Sarría Valencia, Junio 2016 Aplicación de datos LiDAR a la gestión forestal AGRADECIMIENTOS En este apartado, me gustaría agradecer a mi tutor Alfonso Fernández Sarría el tiempo que ha dedicado cada vez que he acudido a su despacho y en cada correo contestado, cuando me surgían las innumerables dudas. En definitiva, sin él, este proyecto no sería posible. También me gustaría agradecer a todos mis compañeros y profesorado, que he tenido el placer de conocer durante el Grado en Ingeniería Geomática y Topografía. Nunca olvidaré ni el apoyo recibido en los momentos duros, ni las risas y las fiestas en los buenos momentos. Para finalizar, me gustaría dar las gracias a mis padres por el cariño recibido durante mi vida. Espero y deseo que sepan lo orgulloso estoy y estaré siempre de ellos. Muchas gracias a todos. Valencia, 14 de junio de 2016 RESUMEN En el ámbito de la gestión forestal, este trabajo tiene como objetivo crear una metodología de procesado e integración de datos LIDAR con imágenes multiespectrales para aplicaciones forestales que requieren una segmentación en rodales. La zona de estudio se encuentra situada en el Parque Natural de la Albufera de Valencia, concretamente entre la Devesa del Saler y la Gola de Pujol. Los datos de partida utilizados son imágenes multiespectrales del satélite Quickbird y datos LIDAR descargados de la página web del CNIG. Se utiliza software libre para los procesos correspondientes a la nube de puntos (FUSION y LASTOOLS) y a la segmentación de los rodales (InterIMAGE). Por otra parte, se utiliza el software ENVI para las operaciones relacionadas con las imágenes multiespectrales y el software ArcGIS como software de apoyo y gestión de la información. En la primera fase del trabajo, se realizan procesos con nubes de puntos tales como recortes, uniones y filtrados hasta conseguir el modelo digital del terreno (MDT) y el modelo digital de superficie (MDS). Una vez obtenidos ambos, podemos calcular el modelo digital de vegetación (MDV o MDSn) donde se representan las alturas de todos los objetos existentes. En la segunda fase del trabajo, se realizan procesos con imágenes multiespectrales tales como recortes, cambios de tamaño de píxel, fusión de imágenes por diferentes métodos, análisis de resultados… Tras escoger el método de fusión con mayor correlación (imagen que combine la resolución espacial de la imagen pancromática y la información espectral de las bandas de la imagen multiespectral), es decir, el método de componentes principales, podemos calcular el Índice de Vegetación de Diferencia Normalizada (NDVI). Este índice representa información sobre el vigor vegetativo y sobre el contenido de humedad de la vegetación. En nuestro caso los valores oscilarán entre -0.011 y 0.540. Valores que cumplen la condición de estar entre -1 y1. En la tercera fase del trabajo, se realiza sobre el MDV distintas segmentaciones mediante el software InterIMAGE y utilizando el algoritmo de Baatz & Shape. La última segmentación se realiza a un archivo resultado de la combinación entre la imagen MDV y la NDVI, cuya solución es la más precisa de todas al utilizar el índice de vegetación. Como se indica en las conclusiones finales, esta última segmentación mejora la identificación de las diferentes cubiertas existentes. Para finalizar, evaluaremos la segmentación realizada mediante la herramienta FUSION CloudMetrics y mediante fotointerpretación. Los diferentes resultados indican que la evaluación ha sido correcta ya que cumplen con los objetivos marcados desde el principio. ABSTRACT In the field of forest management, this research has the aim of create a methodology processing and data integration with LIDAR points and multispectral images for forestry applications requiring polygon segmentation. The study area is located in the Albufera Natural Park of Valencia, specifically between the Devesa del Saler and the Gola de Pujol. The input data used are multispectral Quickbird satellite images and LIDAR data downloaded from the website of CNIG. Free software will be used for processes corresponding to the point cloud (FUSION and LASTOOLS) and polygon segmentation (INTERIMAGE). Moreover, the ENVI software will be used for operations related to the multispectral images and ArcGIS software as software support and management of the information. In the first phase of work, there are processes related with point clouds such as trimmings, unions and filtered until the digital terrain model (DTM) and digital surface model (DSM) are performed. Once we obtained both digital models, we can calculate the model digital model of vegetation (CHM or DSMn) where the heights of all existing objects are represented. In the second phase of work, there are multispectral imaging processes such as trimmings, pixel resizes, image fusion by different methods, analysis of results... After choosing the method of fusion with higher correlation (image that combines the spatial resolution of the panchromatic image and spectral information of multispectral image bands), in this case Principal Components method, we can calculate the Normalized Differential Vegetation Index (NDVI). This index represents information on the vegetative vigor. In our case, these values are right because they are inside -1 and 1 range In the third phase of work, we will be realized different segmentations to CHM image with INTERIMAGE software and using the Baatz & Shape algorithm. The final segmentation will be performed to a file result of the combination between CHM and DMV image, which solution will be the most accurate. As indicated in final conclusions, the last segmentation will improve the identification of different existing decks. Finally, we will evaluate this segmentation by FUSION tool CloudMetrics and by photointerpretation. Different results indicate that the evaluation was successful because they meet the objectives mentioned at the beginning. ÍNDICE 1. INTRODUCCIÓN ____________________________________________________ 1 1.1. OBJETIVOS __________________________________________________________ 3 1.2. LOCALIZACIÓN _______________________________________________________ 4 1.3. ESQUEMA DEL TRABAJO _______________________________________________ 7 1.4. SOFTWARE EMPLEADO ________________________________________________ 8 2. PROCESO CON DATOS LIDAR ________________________________________ 11 2.1. MERGEDATA ________________________________________________________ 11 2.2. POLYCLIPDATA ______________________________________________________ 12 2.3. FILTERDATA ________________________________________________________ 14 2.4. CLIPDATA __________________________________________________________ 15 2.5. CANOPYMODEL _____________________________________________________ 17 2.6. GROUNDFILTER _____________________________________________________ 17 2.7. GRIDSURFACE CREATE ________________________________________________ 18 2.8. DTM2ASCII _________________________________________________________ 19 2.9. CÁLCULO DEL MDSn O MDV ___________________________________________ 20 3. ORTOFOTOGRAFÍA PNOA ___________________________________________ 21 4. PROCESO DE ANÁLISIS DE IMÁGENES MULTIESPECTRALES ________________ 23 4.1. TAMAÑO DE LAS IMÁGENES ___________________________________________ 23 4.2. MÉTODOS DE FUSIÓN ________________________________________________ 25 4.3. CÁLCULO NDVI ______________________________________________________ 33 5. INTRODUCCIÓN AL MÉTODO DE SEGMENTACIÓN _______________________ 37 5.1. ALGORITMO BAATZ-SCHÄPE ___________________________________________ 38 5.2. TRABAJOS DE SEGMENTACIÓN Y CLASIFICACIÓN A REALIZAR _________________ 40 5.3. RECORTE ZONA DE ESTUDIO ___________________________________________ 41 6. RODALIZACIÓN Y EXTRACCIÓN DE LA ALTURA MEDIA DEL RODAL __________ 43 6.1. CARGA DE DATOS DE PARTIDA__________________________________________ 43 6.2. REGLAS DE SEGMENTACIÓN ___________________________________________ 44 6.3. CÁLCULO DE LA ALTURA MEDIA DE CADA SEGMENTO _______________________ 47 6.4. CLASIFICACIÓN DE LOS RODALES ________________________________________ 48 6.5. CARGA DE DATOS EN ARCMAP _________________________________________ 55 6.6. EVALUACIÓN DE LOS RESULTADO OBTENIDOS _____________________________ 56 6.7. ANÁLISIS DE RESULTADOS DE LA MATRIZ DE CONFUSIÓN ____________________ 59 7. TIPIFICACIÓN DE ESTRATOS DE LAS CUBIERTAS _________________________ 61 7.1. REGLAS DE SEGMENTACIÓN ___________________________________________ 61 7.2. ANÁLISIS DE RESULTADOS DE LA MATRIZ DE CONFUSIÓN ____________________ 65 8. DELIMITACIÓN DE RODALES A PARTIR DE LA INTEGRACIÓN DE DATOS LIDAR Y DE DATOS MULTIESPECTRALES Y OBTENCIÓN DE PARÁMETROS ESTRUCTURALES A NIVEL DE RODALES _________________________________________________________ 67 8.1. GEORREFERENCIACIÓN _______________________________________________ 67 8.2. ANÁLISIS DEL NDVI ___________________________________________________ 70 8.3. COMPOSICIÓN DE BANDAS ____________________________________________ 73 8.4. REGLAS DE SEGMENTACIÓN ___________________________________________ 74 8.5. ANÁLISIS DE RESULTADOS DE LA MATRIZ DE CONFUSIÓN ____________________ 78 8.6. EVALUACIÓN DE LOS RESULTADOS MEDIANTE FOTOINTERPRETACIÓN __________ 79 9. DETERMINACIÓN DE PARÁMETROS ESTRUCTURALES DE LOS RODALES ______ 81 10. CONCLUSIONES ___________________________________________________ 95 11. BIBLIOGRAFÍA ____________________________________________________ 97 12. ANEXOS _________________________________________________________ 99 ÍNDICE DE FIGURAS Figura 1. Unidades de suelo en la Albufera. _________________________________________________ 4 Figura 2. Situación zona de estudio 1. ______________________________________________________ 5 Figura 3. Situación zona de estudio 2. ______________________________________________________ 6 Figura 4. Situación zona de estudio 3. ______________________________________________________ 6 Figura 5. Esquema del trabajo. ___________________________________________________________ 7 Figura 6. Software empleado. ____________________________________________________________ 7 Figura 7.Proceso descompresión de .las con LASTOOLS. ______________________________________ 11 Figura 8. Ficheros agrupados en un .txt. ___________________________________________________ 11 Figura 9. Código MergeData. ___________________________________________________________ 12 Figura 10. Resultado de la operación MergeData. ___________________________________________ 12 Figura 11. Código PolyClipData. _________________________________________________________ 12 Figura 12. Límite zona de estudio. ________________________________________________________ 13 Figura 13. Zona de estudio recortada. ____________________________________________________ 14 Figura 14. Outliers. ____________________________________________________________________ 14 Figura 15. Código FilterData. ____________________________________________________________ 15 Figura 16. Código ClipData. _____________________________________________________________ 15 Figura 17. Zona donde se encuentra el outlier. ______________________________________________ 16 Figura 18. Visor LDV donde se observa el outlier. ____________________________________________ 16 Figura 19. Código CanopyModel. ________________________________________________________ 17 Figura 20. Iteraciones del comando GroundFilter. ___________________________________________ 18 Figura 21. Comando GridSurfaceCreate. ___________________________________________________ 18 Figura 22. Comando DTM2ASCII. ________________________________________________________ 19 Figura 23.Resumen de los procesos realizados con FUSION. ___________________________________ 19 Figura 24. Visualización en ArcMap del MDS y MDT. _________________________________________ 20 Figura 25.Visualización en ArcMap del MDSn o MDV. ________________________________________ 20 Figura 26. Ortofotografía descargada del PNOA. ____________________________________________ 21 Figura 27. Tabla estadísticas radiométricas. _______________________________________________ 23 Figura 28. Cambio del tamaño del píxel. ___________________________________________________ 24 Figura 29. Imagen Pancromática y Multiespectral. __________________________________________ 24 Figura 30. Esquema Método de Brovey ___________________________________________________ 25 Figura 31. Imagen Método de Brovey. ____________________________________________________ 25 Figura 32. Espacios de color RGB y HSI. ___________________________________________________ 26 Figura 33.Esquema Método HSI. _________________________________________________________ 26 Figura 34. Imagen Método HSI. __________________________________________________________ 27 Figura 35. Método de Componentes Principales. ____________________________________________ 27 Figura 36. Autovalores de los componentes principales. ______________________________________ 28 Figura 37. Estadísticas de los componentes principales y de la imagen pancromática. ______________ 28 Figura 38. Fórmula Componentes Principales. ______________________________________________ 29 Figura 39. Elementos de la fórmula de Componentes Principales. ______________________________ 29 Figura 40. Estadísticas de los componentes principales y de la imagen pancromática. ______________ 29 Figura 41. Imagen Método Componentes Principales. ________________________________________ 30 Figura 42. Resultado espectral de los métodos utilizados. _____________________________________ 31 Figura 43. Mosaico comparativo de métodos. ______________________________________________ 32 Figura 44. Cambio del tamaño del píxel. ___________________________________________________ 33 Figura 45. Zona de estudio en color verdadero. _____________________________________________ 33 Figura 46. Fórmula para calcular el NDVI. _________________________________________________ 34 Figura 47. Estadísticas NDVI. ____________________________________________________________ 35 Figura 48. Imagen NDVI. _______________________________________________________________ 35 Figura 49. Comparativa NDVI recortado y Ortofoto recortada. _________________________________ 36 Figura 50. Crecimiento de regiones. ______________________________________________________ 37 Figura 51. Factor de fusión. _____________________________________________________________ 38 Introducción 4 1.2. LOCALIZACIÓN El área de trabajo está ubicada entre la zona de la Devesa del Saler y la Gola de Pujol. Estas dos zonas se encuentran localizadas al Este del Parque Natural de la Albufera de Valencia. Tiene una longitud de 5 km y una anchura de 1 km, es decir, presenta una superficie aproximada de 5 km cuadrados. La Devesa de la Albufera se fue formando por los materiales detríticos depositados por la corriente marina. Los materiales de partida fueron también aportados por los ríos y barrancos de la zona, que iniciaron el desarrollo de una barra submarina que posteriormente emergería en forma de flecha, para, finalmente, constituir el cordón dunar que poco a poco fue cerrando el golfo marino existente. Como consecuencia se constituyó el Lago de la Albufera y la Devesa. Por ello, el Lago de la Albufera corresponde al tipo de penilago engendrado por un cordón litoral que se encuentra fijado por la vegetación. En principio, sus aguas fueron salinas pasando paulatinamente a dulces por los aportes fluviales. La presencia de litologías triásicas y cretácicas costeras sirvió como punto de apoyo para las formaciones dunares, así como también del cordón litoral antes mencionado. La presencia de la duna fósil del Perellonet, datada por Sanjaume (1980) en el periodo Tirreniense, indica la antigüedad de la formación. A continuación se muestra una imagen que representa las unidades de suelo existentes en zona de la Albufera: Figura 1. Unidades de suelo en la Albufera. Introducción 5 Se caracteriza, principalmente, por contener grandes dunas ubicadas en la playa que originan un ecosistema donde aparecen tipos de vegetación (pinos, pastos y matorrales generalmente) y fauna (generalmente aves endémicas). Si analizamos la vegetación existente más profundamente, podemos decir que es actualmente pinar de halepensis con algunas manchas de Pinus pinea L. y Pinus pinaster Aiton, con su matorral asociado. El Pinus halapensis L. ha desplazado al Juniperus macrocarpa L y se ha convertido en el dominante de la comunidad. En las dunas interiores, bien cubiertas de vegetación, además de las plantas de la zona de transición, pueden encontrarse plantas trepadoras de hasta 2 m. Finalmente, el viento se constituye como un factor clave en la génesis y dinámica de esta zona, contribuyendo de manera apreciable a la morfología de la Devesa. En el límite Sur de la zona, se encuentra la Gola de Pujol ya mencionada. Se trata del canal más moderno que conecta directamente el lago con el mar. Como dato de interés, cabe destacar, que en 1965 se llevó a cabo un proceso de urbanización que transformó considerablemente el paisaje. En dicho proceso, se destruyeron gran parte de las dunas para poder construir viviendas, vías de comunicación (CV-500), zonas deportivas (campo de golf),… En la actualidad, no se podría repetir este suceso ya que en 1979 se desarrollaron leyes y proyectos para recuperar y regenerar esta zona costera gracias, en gran medida, a la participación de las poblaciones de alrededor. Figura 2. Situación zona de estudio 1. Introducción 6 Figura 3. Situación zona de estudio 2. Figura 4. Situación zona de estudio 3. Introducción 7 1.3. ESQUEMA DEL TRABAJO Figura 5. Esquema del trabajo. LASTOOLS ENVI INTERIMAGE ARCGIS FUSION Figura 6. Software empleado. El proceso comienza con la descompresión de los datos mediante LASTOOLS. Una vez tenemos los datos en .las, realizamos las operaciones pertinentes a la nube de puntos con la ayuda de la ortofoto previamente descargada. Cuando tenemos la nube de puntos preparada, calculamos los MDS y MDT, para posteriormente con ArcGIS realizar la resta y obtener el MDV. Por otra parte, tenemos las imágenes multiespectrales. Con ellas, realizamos diferente métodos de fusión hasta encontrar el más óptimo. Una vez conseguido se calcula el NDVI. Éste índice, junto con el MDV anteriormente citado y la ortofoto serán nuestros datos a la hora de realizar las segmentaciones. Una vez realizadas las segmentaciones que hemos considerado oportunas (rodalización, clasificación de los estratos de masa y la que integra el MDV y el NDVI), las evaluamos mediante matrices de confusión. Para finalizar, calculamos una serie de parámetros estadísticos mediante la herramienta CloudMetrics y realizamos una evaluación visual construyendo planos con cada una de los anteriores parámetros. Introducción 8 1.4. SOFTWARE EMPLEADO LASTOOLS LASTOOLS es un producto de software libre que convierte rápidamente archivos LAS voluminosos en archivos compactos LAZ y viceversa, sin pérdida de información. La compresión realizada por LASTOOLS suele ser más pequeña, y muchas veces más rápida que los compresores genéricos como bz2, gzip, y rar, porque sabe lo que los diferentes bytes en un archivo LAS representan. Otra ventaja de LASTOOLS es que le permite tratar los archivos comprimidos LAZ igual que los archivos estándar de LAS. Puede cargarlos directamente de forma comprimida en su aplicación sin necesidad de descomprimirlos en el disco duro primero. LASTOOLS, fue ganador del premio a la Innovación de Tecnología Geoespacial del Foro Mundial del 2012 en procesamiento LIDAR y subcampeón de producto innovador en INTERGEO 2012 como estándar de facto para la compresión LIDAR. ENVI Consiste en un software de procesamiento y análisis de imágenes desarrollado por EXELIS. ENVI combina procesamiento de las imágenes espectrales con la tecnología de análisis de imágenes mediante una interfaz fácil de utilizar para ayudar a obtener información significativa de las imágenes. En nuestro caso se realizaran tareas como:  Registrar dos o más imágenes.  Fusión de imágenes, máscaras, generación de mosaicos.  Calcular índices de vegetación. FUSION Se trata de un software libre creado por el Servicio Forestal de EEUU para el tratamiento y gestión de datos LIDAR. El sistema de análisis y visualización se compone de dos programas principales, FUSION y LDV (visor de datos LIDAR), y una colección comandos para tareas específicas. La interfaz principal, proporcionada por FUSION, consiste en una ventana de visualización gráfica y una ventana de control. La pantalla presenta todos los datos del proyecto utilizando una pantalla 2D típico de los sistemas de información geográfica. Es compatible con una variedad de tipos de datos y formatos, incluyendo archivos de formas, imágenes, modelos digitales del terreno, modelos de superficie y los datos de retorno LIDAR. LDV proporciona el entorno de visualización 3D para el examen y la medición de subconjuntos de datos espacialmente explícitos. En nuestro caso utilizaremos FUSION para realizar operaciones como uniones, recortes, filtrado de outliers, generación de modelos de superficies, conversión de datos… Introducción 9 ARCGIS Consiste en un software desarrollado por ESRI y que contiene diversas aplicaciones SIG. Las dos aplicaciones principales para nuestro trabajo son ArcMap y ArcCatalog. Cada aplicación cuenta con funciones únicas que se ajustan a las necesidades del usuario. En este caso se trabajará con la versión 10.1 y se realizarán tareas como:  Operaciones de datos ráster para el cálculo del MDV.  Unión de tabla de atributos.  Análisis estadísticos.  Unión de imágenes en una sola.  Visualización de resultados con su posterior edición. INTERIMAGE Se trata de un software de código abierto escrito en C ++ que forma parte de un proyecto de cooperación científica internacional dirigido por el Laboratorio de Visión por Ordenador del Departamento de Ingeniería Eléctrica de la Universidad Católica de Río de Janeiro (PUC-Rio) y por el Instituto Nacional de Investigación Espacial (INPE). Consiste en un software para la clasificación de imágenes usando el principio de segmentación. El concepto es bastante interesante, ya que uno puede modelar su conocimiento de una zona mediante el uso de semántica. Es decir si poseo una imagen satelital y tengo información de la zona, entonces puedo definir clases que me permiten integrar mi conocimiento en la clasificación de las imágenes. En la mayoría de los softwares, el proceso de clasificación se encuentra generalmente determinado por el nivel de reflexión, aunque aquí se usan también estos principios, la novedad de INTERIMAGE es que permite integrar el conocimiento del experto. La única desventaja de InterIMAGE es que solo se pueden usar imágenes y, en nuestro caso, cuando utilizamos imágenes que contienen más de una banda pueden ralentizar considerablemente el proceso. 10 Proceso con datos LiDAR 11 2. PROCESO CON DATOS LIDAR En primer lugar, se han descargado los datos LiDAR comprimidos (.laz) de 2x2 km de la página web del Centro de Descargas del CNIG ( http://centrodedescargas.cnig.es/CentroDescargas ). Una vez obtenidos, con la aplicación laszip.exe descomprimimos los archivos .laz y los convertimos a .las para poder trabajar con FUSIÓN. Figura 7.Proceso descompresión de .las con LASTOOLS. 2.1. MERGEDATA La primera operación a realizar una combinación de varios archivos de nubes de puntos en un solo archivo. Esto se realiza escribiendo la ruta de los ficheros en un único bloc de notas y aplicándoles la operación MergeData para convertir el archivo .txt en un archivo .las ya con el conjunto de datos incorporado. A continuación se muestran las operaciones realizadas en este caso: Figura 8. Ficheros agrupados en un .txt. Proceso con datos LiDAR 12 Figura 9. Código MergeData. El resultado gráfico visualizado mediante la aplicación PDQ sería el siguiente: Figura 10. Resultado de la operación MergeData. 2.2. POLYCLIPDATA Como se puede observar, la zona visualizada es demasiado extensa (10 x 4 km) ya que, en este caso, solo vamos a realizar el estudio de una zona de 5 km de largo y 1km de ancho. Por lo tanto, la segunda operación será realizar un recorte de datos de puntos usando polígonos almacenados en ficheros mediante el comando PolyClipData. La operación a realizar será la siguiente: Figura 11. Código PolyClipData. Proceso con datos LiDAR 13 Como se puede observar en la segunda línea de código, la zona que nos interesa extraer será un shape y ha sido obtenida, previamente, mediante el software ArcGIS. Con el comando Create New Shapefile, obtenemos el límite de la nuestra zona de estudio para agilizar el procesado de datos. Le asignamos el sistema de coordenadas siguiente: ETRS 1989 UTM HUSO 30. Con la ayuda de la ortofoto correspondiente, dibujamos el perímetro de la zona con el comando Editor. El resultado es el que aparece a continuación con la ortofoto incluida: Figura 12. Límite zona de estudio. La tercera línea de código muestra el nombre del resultado del recorte, mientras que la cuarta línea muestra el nombre del fichero al que se le va a realizar la operación del recorte. Proceso con datos LiDAR 20 2.9. CÁLCULO DEL MDSn O MDV Llegado a este punto, se ha finalizado todo el trabajo correspondiente con el programa FUSIÓN. A continuación, seguimos con el software ArcGis. Con el comando Conversion Tools  ASCII to Raster convertimos los modelos a formatos ráster. El resultado es el siguiente: Figura 24. Visualización en ArcMap del MDS y MDT. Una vez obtenidos los modelos anteriores, necesitamos saber las alturas de los objetos ubicados en nuestra zona de estudio. Para ello, se calcula el Modelo Digital de Superficies Normalizado (MDSn) también llamado Modelo Digital de Vegetación (MDV). Utilizaremos la herramienta Raster Calculator para conseguirlo. Este proceso se basa en la sencilla operación de restarle al MDS el MDT. El resultado es el siguiente: Figura 25.Visualización en ArcMap del MDSn o MDV. Ortofotografía PNOA 21 3. ORTOFOTOGRAFÍA PNOA En 2004, se crea el Plan Nacional de Ortofotografía Aérea (PNOA). Consiste en un proyecto cofinanciado y cooperativo entre la Administración General del Estado (AGE) y las comunidades autónomas que se enmarca dentro del Plan Nacional de Observación del Territorio (PNOT), siendo coordinado por el Instituto Geográfico Nacional (IGN) y el Centro Nacional de Información Geográfica (CNIG). Tiene como objetivo la obtención de ortofotografías digitales para todo el territorio nacional, incluyendo: el vuelo fotogramétrico, apoyo de campo, aerotriangulación y el modelo digital de elevaciones. El proyecto se encuentra en continua evolución, adaptándose a las necesidades de los usuarios y al desarrollo de nuevas tecnologías. En 2009, se planteó la obtención de Modelos Digitales de alta precisión, obtenidos por tecnología LIDAR, para la realización de cartografía de áreas de inundación, proyectos de carreteras, inventarios forestales, etc. Se han descargado la imagen desde el Centro de Descargas del Centro de Información Geográfica (CNIG) en formato .ecw que nos servirá como cartografía base. Figura 26. Ortofotografía descargada del PNOA. 22 Proceso de análisis de imágenes multiespectrales 23 4. PROCESO DE ANÁLISIS DE IMÁGENES MULTIESPECTRALES 4.1. TAMAÑO DE LAS IMÁGENES Partimos con dos imágenes QuickBird-2 (una de ellas pancromática y otra multiespectral tomadas el 16 de noviembre de 2015) que tenemos intención de fusionar para poder mejorar propiedades tales como, la resolución espacial de una imagen multiespectral, mejorar los procesos de segmentación (tarea que también realizaremos en este trabajo) y mejorar la clasificación de las texturas en imágenes. Las dos imágenes, aparentemente, representan la misma zona pero son de diferente tamaño (distintas filas y columnas) y resolución espacial. La imagen pancromática tiene una superficie de 7.915 km2 (3944 columnas y 5575 filas con un tamaño de píxel de 0.6 metros. La imagen multiespectral tiene una superficie de 7.403 km2 (3944 columnas y 5575 filas con un tamaño de píxel de 2.4 metros. En la siguiente tabla se resumen los parámetros radiométricos: Imagen Banda Resolución espectral (µm) Pancromática Pan 0.45 - 0.90 Multiespectral Azul 0.45 - 0.52 Verde 0.52 - 0.60 Rojo 0.63 - 0.69 Infrarrojo 0.76 - 0.90 Figura 27. Tabla estadísticas radiométricas. La primera tarea a realizar por tanto será la del cambio de resolución espacial de la imagen multiespectral de 2.4 m a 0.6 m que es el mismo que el de la imagen pancromática. Proceso de análisis de imágenes multiespectrales 24 Figura 28. Cambio del tamaño del píxel. La segunda tarea será cambiar el tamaño de la imagen pancromática con el comando Resize Data por el tamaño de la imagen multiespectral con la resolución modificada anteriormente. Observamos el nuevo tamaño y vemos que tiene 3907 filas y 5136 columnas. Sabiendo esto podremos comprobar al realizar el siguiente cambio de tamaño si lo hemos realizado bien o no, ya que las dos imágenes tienen que tener el mismo número de filas y columnas para poder realizar la fusión. La tercera tarea será realizar otro cambio de tamaño, pero esta vez de la imagen multiespectral por el tamaño de la imagen pancromática recortada anteriormente. El resultado de las dos imágenes se muestra a continuación: Figura 29. Imagen Pancromática y Multiespectral. Ambas imágenes tiene el mismo número de filas y columnas y la misma resolución espacial. Proceso de análisis de imágenes multiespectrales 25 4.2. MÉTODOS DE FUSIÓN Ya estamos en disposición de aplicar los respectivos métodos de fusión. Para esta tarea se aplicarán un conjunto de métodos y posteriormente se seleccionará el más óptimo. 1. Método de Brovey Este método sigue el siguiente esquema: Figura 30. Esquema Método de Brovey . Obtenemos la imagen fusionada realizando el producto entre cada banda multiespectral y la banda pancromática y dividiendo por la suma de las bandas de la imagen multiespectral. Con las herramientas File /Save File as / Envi Standard agrupamos las tres imágenes en un fichero. También podemos realizar la fusión utilizando un método directo de Envi con la herramienta Transform\Image Sharpening\Color Normalized (Brovey). El resultado de dicha fusión es el siguiente: Figura 31. Imagen Método de Brovey. Proceso de análisis de imágenes multiespectrales 26 2. Método HSI Para visualizar una imagen digital en color se necesita un espacio de representación de los colores compatible con los dispositivos de visualización utilizados. Gracias al espacio de coordenadas RGB, donde las coordenadas representan la proporción de cada uno de los colores primarios (rojo, verde y azul, respectivamente) se consigue la representación de cada punto de la imagen en color. Sin embargo, según el punto de vista humano, se debería representar dichas coordenadas en función de propiedades, tales como: tono saturación e intensidad. Este se puede conseguir con el espacio HSI (Hue, Saturation and Intensity). En la siguiente imagen se muestran los dos espacios de coordenadas mencionados anteriormente: Figura 32. Espacios de color RGB y HSI. Este método sigue el siguiente esquema: Figura 33.Esquema Método HSI. Proceso de análisis de imágenes multiespectrales 27 En este método, se realiza la transformación de tres bandas del espacio de color RGB al espacio de color HSI. Se sustituye el componente Intensidad (Lightness) por la imagen pancromática y se realiza la transformación inversa del espacio HLS a RGB con la herramienta Transform /Image Sharpening/HSV. El resultado de dicha fusión es el siguiente: Figura 34. Imagen Método HSI. 3. Método de Componentes Principales Se trata de un método que genera nuevas bandas que representan la mayor cantidad de variabilidad de los datos utilizando el menor número de componentes posible. El primer componente contiene más información que el segundo, el segundo más que el tercero, y así sucesivamente. El primer componente suele asociarse a la intensidad de la imagen, por lo tanto, cuando sea sustituido por la imagen pancromática (mayor resolución) no aparecerán diferencias espectrales considerables. Dicho método cumplirá el siguiente esquema: Figura 35. Método de Componentes Principales. Proceso de análisis de imágenes multiespectrales 28 Comenzamos ejecutando el comando Transform / Principal Components / Forward PC Rotation / Compute New Static and Rotate. Una vez ejecutado el comando, aparece una ventana que contiene los autovalores de los componentes principales como se puede observar en la siguiente imagen: Figura 36. Autovalores de los componentes principales. Como se puede apreciar en la imagen anterior, la mayoría de la información se encuentra almacenada en el primer componente principal. Sabiendo esto, se hallan las estadísticas (media y desviación típica) de la imagen pancromática y del primer componente principal. El resultado se muestra en la siguiente página. Figura 37. Estadísticas de los componentes principales y de la imagen pancromática. Proceso de análisis de imágenes multiespectrales 29 MEDIA DESVIACIÓN PAN 255.626953 190.129919 CP1 0 358.907810 Para poder sustituir las imágenes, se debe ajustar primero la imagen pancromática respecto a los datos del CP1. Siguiendo la siguiente formulación: Panajustada = ( Pan * a ) +b Figura 38. Fórmula Componentes Principales. Siendo: a = σcp1 / σPan = 1.887697 b = µCP1 – ( µPan * a ) = -482.546232 Figura 39. Elementos de la fórmula de Componentes Principales. Con la herramienta Band Math ((b1 *a) + b, siendo b1 = Pan) y obtenemos las estadísticas que se muestran a continuación. Figura 40. Estadísticas de los componentes principales y de la imagen pancromática. Proceso de análisis de imágenes multiespectrales 36 Como posteriormente se verá en el apartado de segmentación, la zona de estudio es demasiado extensa, de manera que procedemos a realizar un recorte de una zona que contenga varios tipos de cubiertas, utilizando ArcMap y su herramienta Spatial Analyst / Extraction / Extract By Mask. La máscara será la imagen MDV recortada. Lo hacemos asía puesto que nos interesa tener tanto el MDV como el NDVI con las mismas medidas a la hora de realizar la segmentación. Como se verá posteriormente, estas medidas no deben sobrepasar los 500 metros cuadrados para que la segmentación no sea errónea. A continuación, se muestra el resultado del recorte y su correspondiente imagen comparativa de la ortofoto recortada. Figura 49. Comparativa NDVI recortado y Ortofoto recortada. Introducción al método de segmentación 37 5. INTRODUCCIÓN AL MÉTODO DE SEGMENTACIÓN La segmentación de imágenes divide la imagen en sus partes constituyentes hasta un nivel de división en el que se aíslen las regiones u objetos de interés. El objetivo de nuestra segmentación será el de realizar una rodalización de las distintas cubiertas existentes en la zona de estudio. Los algoritmos empleados en este proceso pueden basarse en dos factores: en la discontinuidad o en la similitud entre los niveles de gris de los píxeles vecinos. Si hablamos de discontinuidad, el objetivo es dividir la imagen teniendo en cuenta los cambios bruscos de niveles de gris. Por ejemplo, realizando una detección de puntos aislados, de líneas o de bordes. Si hablamos de similitud, el objetivo es dividir la imagen teniendo en cuenta las zonas que contengan valores similares. Existen varios procesos como el crecimiento de regiones y la umbralización. Figura 50. Crecimiento de regiones. El algoritmo que vamos a emplear para nuestro trabajo se basa en el crecimiento de regiones. Este método consiste en la suposición de propiedades similares entre píxeles adyacentes. Se desarrolla en 4 pasos: 1. Define un grupo de puntos previos, también llamados “semillas”. 2. De forma iterativa, se analizan los vecinos de cada “semilla”. 3. Se pueden incorporar nuevos píxeles vecinos a la región si cumplen alguna condición propuesta por el usuario. 4. El proceso finaliza cuando no se produzcan más cambios o cuando se llegue a un número máximo de iteraciones. El inconveniente que presentan, los algoritmos basados en crecimientos de regiones, es la gran carga de procesamiento que tiene que soportar el ordenador utilizado debido a que se trata de un proceso iterativo. Introducción al método de segmentación 38 5.1. ALGORITMO BAATZ-SCHÄPE El algoritmo empleado en nuestro caso será el Baatz-Schäpe. Consiste en un proceso iterativo, que busca minimizar la heterogeneidad promedio de los objetos de una imagen resultantes. Inicialmente, todos los píxeles de la imagen se consideran como segmentos iniciales y cada iteración calcula el aumento de la diversidad, también llamado factor de fusión, como resultado de una posible fusión existente entre cada par de segmentos adyacentes. El algoritmo toma una decisión para elegir la fusión de uno o de varios pares de segmentos que cumplen con los criterios de heterogeneidad. Existen cuatro alternativas que se presentan a continuación:  Decisión "vecino arbitrario". Dos segmentos adyacentes se funden tan pronto como se compruebe que cumplen con los criterios de heterogeneidad. Esta selección de pares de segmentos agregados está fuertemente influenciada por el orden en que los segmentos son procesados por el algoritmo.  Decisión "mejor ajuste" (Best Fitting). Establece que un segmento debería fusionarse con el segmento adyacente por de cumplir con los resultados del criterio de heterogeneidad y por tener un bajo incremento de la heterogeneidad entre todos sus vecinos. La secuencia de visitas en este caso tiene un menor impacto sobre el resultado.  Decisión "mejor ajuste mutuo" (Mutual Best Fitting). La fusión se lleva a cabo sólo si la relación de mayor similitud es mutua. Es decir, dado un segmento A, busca un segmento vecino B que mejor cumple los criterios de homogeneidad. Luego busca para B el objeto vecino C, para el que B cumple mejor los criterios de homogeneidad. Si C = A entonces fusiona los objetos, en caso contrario repite el proceso con respecto a B para A y C para B. Esta secuencia tiene un impacto aún menor que el anterior en el resultado final de la segmentación.  Decisión "mejor ajuste global" (Global Mutual Best Fitting). Establece que sólo el par de mejores vecinos mutuos puede dar como resultado un menor aumento de la heterogeneidad entre los pares de vecinos de toda la imagen. En esta decisión el resultado final de la segmentación es independiente de la secuencia de las visitas. A continuación se muestra la formulación referente a este factor: Figura 51. Factor de fusión. Introducción al método de segmentación 39 El factor de fusión, se define por la suma ponderada de un componente relacionado con la heterogeneidad espectral (hcor) y otra relacionada con la heterogeneidad morfológica (hforma). La importancia relativa de los componentes se define por un peso (wcor) y el cálculo de los componentes sobre la base de la diferencia entre el potencial objeto generado (obj3) y la suma independiente del objeto (obj1 y obj2) como se muestra en la segunda ecuación. La heterogeneidad espectral (hcor) viene dada por la suma ponderada de la desviación estándar de los valores de los píxeles que componen el segmento. Un peso está asociado con cada banda espectral para expresar su importancia relativa. La heterogeneidad morfológica (hforma) se compone de dos elementos diferentes: compacidad y suavidad. La compacidad (Cmp), es la relación de la longitud del borde (b) del segmento y la raíz cuadrada de su área (n). Figura 52. Compacidad. La suavidad (Svd) es la relación entre la longitud del borde (b) del segmento y la longitud del borde (bbox) su cuadro delimitador mínimo Figura 53. Suavidad. Un peso establecido por el usuario expresa la importancia relativa de la compacidad y la suavidad en la composición de la heterogeneidad morfológica. Este valor debe ser inferior a un umbral dado, llamado escala, de modo que la agregación se puede lograr. El proceso se repite hasta que ya no es posible llevar a cabo las fusiones. Los parámetros tales como la relevancia de cada banda espectral y la importancia relativa de la forma y del color, y entre compacidad y suavidad se pueden ajustar con el fin de lograr un mejor resultado en la segmentación Cabe señalar que para llevar a cabo un crecimiento de regiones, cada objeto se selecciona sólo una vez cada iteración. Además, la selección de los objetos se lleva a cabo con el fin de seleccionar objetos relativamente distantes uno del otro de acuerdo a la ubicación en la imagen. Introducción al método de segmentación 40 A diferencia de otros algoritmos basados en crecimiento de regiones, el Baatz & Shape tiene en cuenta tanto la respuesta espectral como la morfológica de cada una de las regiones a la hora de fusionarlas. Los atributos morfológicos son importantes porque trabajamos con una imagen (MDV o MDSn) que solo posee una banda (intensidad). Por lo tanto, consideramos obligatorio tener en cuenta las propiedades espaciales también. A pesar de ser un algoritmo más complejo de ajustar que el resto, tiene la gran ventaja de producir soluciones más precisas. 5.2. TRABAJOS DE SEGMENTACIÓN Y CLASIFICACIÓN A REALIZAR En esta apartado se describe brevemente todas las tareas que vamos a llevar a cabo con el software InterIMAGE relacionados con el proceso de segmentación.  Rodalización de la cubierta y extracción de la altura media del rodal. El objetivo de esta tarea es segmentar la imagen hasta que nos permita identificar los diferentes rodales y calcular su altura media. Posteriormente realizaremos una clasificación de los diferentes tipos de cubiertas a partir de los rodales anteriores. Para finalizar, realizaremos una evaluación estadística calculando la matriz de confusión. Para ello, dibujaremos regiones que sepamos con certeza a la clase que pertenecen teniendo como ayuda la ortofotografía del PNOA que tenemos como dato de partida.  Tipificación de estratos de la masa. El objetivo de esta tarea es segmentar la imagen a nivel objeto hasta que nos permita realizar una distinción de las diferentes cubiertas vegetales de nuestra zona de estudio.  Delimitación de rodales a partir de la integración de datos LIDAR y de datos multiespectrales y obtención de parámetros estructurales a nivel de rodales. El objetivo de esta tarea es mejorar las precisiones a la hora de identificar rodales, introduciendo el índice NDVI calculado a partir de las imágenes multiespectrales QuickBird. Utilizaremos un parámetro de escala más pequeño que en los pasos anteriores para que sea menos complicado rechazar aquellos segmentos que no cumplan las condiciones introducidas. Introducción al método de segmentación 41 5.3. RECORTE ZONA DE ESTUDIO El software que vamos a emplear (InterIMAGE) para las tareas de segmentación tiene el inconveniente de producir errores en el procesamiento si se trabaja con zonas extensas. Por lo tanto vamos a proceder a realizar un recorte del MDV como ya se ha comentado en apartados anteriores. Las coordenadas serán las siguientes: Xmax = 731053 m Xmin = 730647 m Ymin = 4360223 m Ymax = 4360678 m Figura 54. MDV recortado. También, vamos a proceder a realizar un recorte a la ortofoto, 100 metros superior por cada coordenada con respecto al MDV. Esto es debido a que puede ser interesante saber que ocurre en los bordes de la imagen MDV. Las coordenadas serán las siguientes: Xmax = 731153 m Xmin = 730547 m Ymin = 4360123 m Ymax = 4360778 m Figura 55. Ortofoto recortada. 42 Rodalización y extracción de la altura media del rodal 43 6. RODALIZACIÓN Y EXTRACCIÓN DE LA ALTURA MEDIA DEL RODAL 6.1. CARGA DE DATOS DE PARTIDA En primer lugar, abrimos InterIMAGE y creamos un nuevo proyecto con el comando File / New Project: Figura 56. Ventana nuevo proyecto en InterIMAGE. Rellenamos los apartados como se indica en la figura anterior. Solo podemos asignarle la opción Default Image a una imagen (en nuestro caso el MDV). Será la única imagen que el programa procesará. Añadimos tanto el MDV como la ortofoto recortada para poder realizar una interpretación visual. Figura 57.Ventana nuevo proyecto en InterIMAGE. Rodalización y extracción de la altura media del rodal 44 6.2. REGLAS DE SEGMENTACIÓN En esta primera fase de segmentación, solo nos interesa realizar una rodalización con nuestros criterios de la escena en general. Por lo tanto, solo nos hace falta un nodo hijo que llamaremos Cover (Cubierta). Para ello, pulsamos con el botón derecho sobre el nodo Scene y ejecutamos el comando Insert Child / Node, dándonos el siguiente resultado: Figura 58.Nodo Scene con su nodo hijo Cover. Ahora seleccionamos el nodo Cover y editamos sus valores en la ventana Node Editor. Le asignamos el color verde y el operador TopDown, citado anteriormente, TA_Baatz_Segmenter. La secuencia del algoritmo será la siguiente: 1. Lee la imagen. 2. Realiza la segmentación. 3. Aplica la regla de decisión implementada por el usuario 4. Genera el archivo de salida. Los parámetros que modificaremos son:  Scale Parameter (sp): valor del atributo de la escala (cualquier valor positivo).  Color Weight: peso del atributo de color de Baatz-Schäpe. (0-1).  Compactness Weight: peso del atributo de compacidad de Baatz-Schäpe. (0-1). El resto de parámetros se dejarán por defecto. Guardamos el proyecto con el comando File / Save Project y lo ejecutamos con el icono . El resultado será una capa vectorial, llamada Result, que representará la segmentación sobre la imagen seleccionada como se indica en la página siguiente. Rodalización y extracción de la altura media del rodal 45 Figura 59. Ejemplo resultado de segmentación. En la ventana Layers, podemos modificar el color de fondo y el de los bordes de los segmentos. Si además queremos visualizar la ortofoto de fondo podemos hacerlo desde la ventana Layers también. Para nuestra zona estudio hemos realizado diferentes pruebas para analizar cual se adapta mejor a las características del terreno: 1ª PRUEBA: Sp = 20 wcolor = 0.5 wcmpct = 0.5 Figura 60. 1ª Prueba de segmentación. Se puede observar que la segmentación no ha sido correcta ya que los polígonos no se adaptan para nada a las distintas cubiertas y objetos existentes en el terreno. Rodalización y extracción de la altura media del rodal 52 DISTANCIA EUCLÍDEA NORMALIZADA ENTRE CLASES EDIFICIOS ÁRBOLES MATORRAL SUELO DESNUDO EDIFICIOS 0,000 11,997 18,363 19,845 ÁRBOLES 11,997 0,000 6,728 3,330 MATORRAL 18,363 6,728 0,000 1,614 SUELO DESNUDO 19,845 8,330 1,614 0,000 Figura 73. Tabla Distancia euclídea normalizada entre clases. Se puede observar como los resultados de las tablas anteriores siguen una lógica coherente en cuanto a las distancias entre clases. Por ejemplo, si un rodal de matorral tiene una altura media de 1,640 metros, la diferencia entre ésta y el suelo sería de 1.614-1.575 metros, es decir, el suelo tendría una altura media de 0.019 metros. Este valor está dentro del rango de la clase suelo como se puede ver en las tablas anteriores, lo que nos permite afirmar que las clases asignadas y sus polígonos han sido correctamente seleccionados. Estos resultados los tendremos en cuenta a continuación para realizar la clasificación. Una vez realizado los cálculos previos a la clasificación, ya podemos seguir trabajando con InterIMAGE. Creamos un nuevo proyecto y construimos una red que contenga al nodo Cover y a sus nodos hijos (las clases que ya hemos comentado antes) como se muestra a continuación: Figura 74. Nodo Cover con sus nodos hijos. Definimos los parámetros y el operador a todos los nodos hijos como se muestra en la página siguiente. Rodalización y extracción de la altura media del rodal 53 Figura 75. Parámetros de los nodos hijos. Volvemos al nodo Cover para configurar las reglas de decisión TopDown Decision Rule. Para introducir las decisiones, seleccionamos la ventana Bottom Up Decision Rule y vamos a emplear la herramienta Selection como se indica a continuación: Figura 76. Herramienta Selección. Rodalización y extracción de la altura media del rodal 54 El cuadro final donde se muestran las alturas asignadas a cada clase quedará de la siguiente forma: Figura 77. Decision Rule completo del nodo Cover. Una vez finalizado, le damos a OK, salvamos el proyecto y ejecutamos. Dándonos el siguiente resultado: Figura 78. Resultado de la rodalización con InterIMAGE. Rodalización y extracción de la altura media del rodal 55 Se puede observar en la imagen anterior, como la clase edificios se representa con color rosa, los árboles en verde oscuro, el matorral en verde claro y el suelo en un tono marrón. Como solo tenemos en cuenta las alturas de las cubiertas, obviamente, aparecen fallos en la clasificación. No es preocupante puesto que en las segmentaciones posteriores se solucionan en gran medida esta clase de errores. Hay que recordar que el objetivo de esta segmentación, no es otro que rodalizar la zona de una forma rápida sin tener muy en cuenta la precisión. 6.5. CARGA DE DATOS EN ARCMAP Una vez finalizado el proceso, nos gustaría poder analizar cada clase por separado y no en conjunto. Por lo tanto, la mejor forma de hacerlo es exportando cada clase en formato shape. Si nos vamos a la ventana Layers, pulsamos en botón , seleccionamos la pestaña Selection y seleccionamos la clase a cargar. Tendremos en cuenta la opción Stage (Bottom Up), el keyname que consideremos y colores de relleno y borde. Para finalizar, pulsamos el botón para cargar la capa. Una vez hemos realizado la carga de todas nuestras capas, seleccionamos las clases una por una y pulsamos el icono de Export as Shapefile . Se abrirá una ventana emergente, donde tendremos que introducir la expresión de la media de las alturas que es la que queremos que aparezca en la tabla de atributos del shapefile. Como se muestra a continuación: Figura 79. Expresión promedio de altura Rodalización y extracción de la altura media del rodal 56 6.6. EVALUACIÓN DE LOS RESULTADO OBTENIDOS Para poder evaluar los resultados, vamos a realizar una comparativa entre los resultados de InterIMAGE y los rodales realizados con ArcMap. Esta comparativa se llevará acabo construyendo una matriz de confusión. Este método consiste en una matriz donde cada columna representa el número de predicciones de cada clase (clases calculadas mediante InterIMAGE), mientras que cada fila representa a las instancias en la clase real (rodales realizados en ArcMap). Uno de los beneficios de las matrices de confusión es que facilitan ver si el sistema está confundiendo dos clases. Si en los datos de entrada, el número de muestras de clases diferentes cambia mucho la tasa de error del clasificador no es representativa de lo bien o mal que realiza la tarea el clasificador. Como tenemos que trabajar con segmentos del mismo tamaño, tendremos que recortar la zona segmentada con los rodales dibujados como muestras evaluatorias. Emplearemos la orden Analysis Tools / Extract / Clip. El resultado de ese recorte, con la ortofoto de fondo, se muestra a continuación: Figura 80.Muestras evaluatorias. . Donde recordemos que la clase 1 es árboles, la clase 2 es matorral, la clase 3 es suelo y la clase 4 edificios Rodalización y extracción de la altura media del rodal 57 Para poder calcular la matriz de confusión, comentada anteriormente, necesitamos saber la superficie de las clases elegidas. Si abrimos la tabla de atributos de la capa resultado, podemos seleccionar las áreas de las clases y la superficie total que será muy importante para el cálculo de la matriz de confusión. Para ello, nos vamos al campo Shape_Area y seleccionamos la opción Statistics. El resultado de la superficie total aparece en el apartado Sum, como se muestra a continuación: Figura 81. Estadísticas de los rodales. Teniendo en cuenta esto último, ya podemos ir rellenando la matriz, dándonos el resultado final siguiente: MATRIZ DE CONFUSIÓN real seg ÁRBOLES MATORRAL SUELO DESNUDO EDIFICIOS TOTAL FIABILIDAD USUARIO ÁRBOLES 1469,980 36,798 1,274 21,183 1529,235 96,13% MATORRAL 19,196 745,358 382,133 0,000 1146,687 65,00% SUELO DESNUDO 0,000 4,811 2381,184 0,000 2385,995 99,80% EDIFICIOS 425,200 123,893 14,266 522,263 1085,622 48,11% TOTAL 1914,376 910,86 2778,857 543,446 6190,847 FIABILIDAD PRODUCTOR 76,79% 81,83% 85,69% 96,10% FIABILIDAD GLOBAL 83,27% KAPPA 0,65 BUENA Figura 82. Matriz de confusión y fiabilidades. Rodalización y extracción de la altura media del rodal 58 Donde:  Fiabilidad del usuario: índice que muestra la probabilidad de que una superficie clasificada dentro de una clase pertenezca realmente a ella. Se calcula de la siguiente forma: 𝑭𝑼=𝑺𝒄𝒍𝒂𝒔𝒆 𝑺𝒕𝒐𝒕𝒂𝒍_𝒓𝒆𝒂𝒍 Figura 83. Fórmula Fiabilidad del usuario .  Fiabilidad del productor: índice que muestra la proporción de una clase que esté correctamente clasificada. Se calcula de la siguiente forma: 𝑭𝑷=𝑺𝒄𝒍𝒂𝒔𝒆 𝑺𝒕𝒐𝒕𝒂𝒍_𝒔𝒆𝒈𝒎𝒆𝒏𝒕𝒂𝒅𝒂 Figura 84. Fórmula Fiabilidad del productor.  Fiabilidad Global: índice que muestra la exactitud y calidad de la matriz de clasificación. Se calcula de la siguiente forma: 𝑭𝑮=∑𝑫𝒊𝒂𝒈𝒐𝒏𝒂𝒍 𝑺𝒕𝒐𝒕𝒂𝒍_𝒄𝒍𝒂𝒔𝒊𝒇𝒊𝒄𝒂𝒅𝒂 Figura 85. Fórmula Fiabilidad Global. En nuestro caso la 𝑆𝑡𝑜𝑡𝑎𝑙_𝑐𝑙𝑎𝑠𝑖𝑓𝑖𝑐𝑎𝑑𝑎 es igual a 6147.539 m2.  Coeficiente Kappa: Este estadístico es una medida de la diferencia entre la exactitud lograda en la clasificación con un clasificador automático (InterIMAGE) y la probabilidad de lograr una clasificación correcta con un clasificador aleatorio. Se calcula, para nuestro caso, de la siguiente forma: 𝒌=(∑𝑫𝒊𝒂𝒈𝒐𝒏𝒂𝒍−∑(𝑺𝒄𝒍𝒂𝒔𝒆_𝒕𝒐𝒕_𝒓𝒆𝒂𝒍 . 𝑺𝒄𝒍𝒂𝒔𝒆_𝒕𝒐𝒕_𝒔𝒆𝒈 𝑺𝒕𝒐𝒕_𝒄𝒍𝒂𝒔 ) 𝒇 𝒊=𝟏 𝒇 𝒊=𝟏 ) (𝑺𝒕𝒐𝒕_𝒄𝒍𝒂𝒔 − ∑(𝑺𝒄𝒍𝒂𝒔𝒆_𝒕𝒐𝒕_𝒓𝒆𝒂𝒍 . 𝑺𝒄𝒍𝒂𝒔𝒆_𝒕𝒐𝒕_𝒔𝒆𝒈 𝑺𝒕𝒐𝒕_𝒄𝒍𝒂𝒔 ) 𝒇 𝒊=𝟏 ) Figura 86. Fórmula Coeficiente Kappa. Rodalización y extracción de la altura media del rodal 59 El resultado de este coeficiente debe de estar en el rango entre 0 y 1. Los rangos de calidad se pueden clasificar de la siguiente manera: COEFICIENTE k SIGNIFICADO < 0,20 POBRE 0,21 - 0,40 DÉBIL 0,41 - 0,60 MODERADA 0,61 - 0,80 BUENA 0,81 - 1 MUY BUENA Figura 87. Tabla valores coeficiente kappa. 6.7. ANÁLISIS DE RESULTADOS DE LA MATRIZ DE CONFUSIÓN Si observamos el resultado del coeficiente kappa nos da un resultado bueno, con una fiabilidad global del 83,27 %. El resultado que más destaca en cuanto a la fiabilidad es la de la clase árboles. Tiene una fiabilidad del usuario muy alta (96,13 %) y la fiabilidad del productor más baja de las cuatro clases (76,79 %), pese a que se trata también de una fiabilidad alta. Esto quiere decir que se ha clasificado como árboles el 76,79 % pese a que el 96,13 % pertenece realmente a esta clase. Por lo tanto, casi el 23% de la clase considerada previamente como árboles ha sido clasificado como edificios. La clase edificios es el caso contrario. Tiene una fiabilidad del usuario relativamente baja (48,11 %), sin embargo una fiabilidad del productor del 96,10 % (la más alta). Esto quiere decir que se ha clasificado como edificios el 96,10 % pese a que el 48,11 % pertenece realmente a esta clase. Por lo tanto, el 51,89 % de la clase considerada previamente como edificios ha sido clasificado como árboles (error por exceso). Este último problema es debido a que los edificios tienen forma de escalón y están rodeados de zonas de árboles cuya altura promediada es muy alta. Esto provoca, que algunos árboles superen las zonas más bajas de estos edificios. De ahí que algunos segmentos de edificios se hayan clasificado como árboles. Es muy complicado realizar clasificaciones teniendo solamente en cuenta la altura promediada, especialmente en los casos de rodales. Por último, la mayor dificultad del proceso es la distinción entre las clases matorral y suelo desnudo. Es muy confusa, incluso teniendo la ortofoto como referencia. En conclusión, el mayor rendimiento lo podemos obtener para la identificación de grupos de árboles para posteriormente cálculos de superficie, biomasa,… 60 Tipificación de estratos de las cubiertas 61 7. TIPIFICACIÓN DE ESTRATOS DE LAS CUBIERTAS Consiste en una tarea similar a la anterior, con la variante de realizar la segmentación a nivel objeto y no a nivel rodal, es decir, el componente morfológico tendrá mayor importancia que el color... Considerando la clasificación de las alturas realizada en el apartado anterior, vamos a aplicar las siguientes reglas para clasificar las siguientes cubiertas:  Edificios (Buildings): alturas superiores a 14 metros.  Arbolado (Tree): alturas comprendidas entre 14 y 3 metros.  Matorral (Grass): alturas comprendidas entre 3 y 0.2 metros.  Suelo (Soil): alturas inferiores a 0.2 metros. 7.1. REGLAS DE SEGMENTACIÓN Creamos un nuevo proyecto y construimos una red que contenga al nodo Region y a sus nodos hijos (las clases que ya hemos comentado antes) como se muestra a continuación: Figura 88. Nodo Region con sus hijos.  Definimos los parámetros y el operador a todos los nodos hijos de la siguiente forma: Figura 89. Parámetros nodos hijos del nodo Region. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 68 Para realizar el proceso de forma correcta, pulsamos la opción Fit to Display dentro de Georreferencing, y posteriormente pulsamos Viewer para poder visualizar ambas imágenes a la vez. Se procede a insertar puntos tanto en el NDVI (imagen a corregir) como en la ortofoto (puntos homólogos) como se muestra a continuación: Figura 99. Puntos de apoyo. Después de insertar los puntos, pulsamos el botón View Link Table para poder controlar los residuos y el error medio cuadrático. Figura 100. Residuos y error medio cuadrático. Consideramos como error máximo 1 metro. Observamos que los residuos entran dentro de lo previsto, por lo tanto se puede guardar el archivo con la geometría modificada mediante el comando Rectify como se muestra en la página siguiente. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 69 Figura 101. Interpolación vecino más cercano. El resultado final es el siguiente: Figura 102. NDVI bien georreferenciado. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 70 Si comprobamos las propiedades de las imágenes MDV y NDVI pulsando el botón derecho sobre ellas observamos que ambas coinciden con los siguientes valores: Figura 103. Propiedades del MDV y NDVI. 8.2. ANÁLISIS DEL NDVI Tras realizar la georreferenciación, tenemos que analizar los valores del NDVI recortado. Analizando bien el índice con ArcMap, vemos que si lo clasificamos en cuatro clases (edificios, árboles, matorral y suelo desnudo) se le asignan los siguientes rangos (max, min): EDIFICIOS ÁRBOLES MATORRAL SUELO NDVI MIN -0,011 0,317 0,227 0,131 NDVI MAX 0,131 0,540 0,317 0,227 NDVI MEDIO 0,009 0,429 0,272 0,179 Figura 104. Valores del NDVI para las diferentes clases. Estos rangos han sido escogidos tras realizar pruebas visuales mediante ArcGIS y modificando los valores de los umbrales de la imagen NDVI hasta poder discriminar claramente cada cubierta. Dentro del rango de edificios, podemos crear otra clase llamada Otros que haga referencia a construcciones tales como aparcamientos, carreteras o piscinas y también a una laguna que se encuentra situada el noroeste de nuestra zona. Estos elementos se pueden apreciar en la siguiente página. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 71 Figura 105. Elementos de la clase "Otros". El rango de la clase Otros sería el siguiente: NDVI MIN NDVI MAX NDVI MEDIO OTROS -0,011 0,121 0,009 Figura 106. Valores NDVI de la clase "Otros". Como se verá en el apartado 8.3, dado que en esta nueva clase no coinciden los parámetros MDV y NDVI a la vez con el resto de clases no habrán problemas a la hora de realizar la segmentación. Con ArcMap, podemos asignarle diferentes colores a los rangos de las cuatro primeras clases, como se muestra en la página siguiente. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 72 Figura 107. NDVI clasificado. Donde: Figura 108. Leyenda del NDVI. Esta información, junto con las alturas calculadas en pasos anteriores, será la que introduciremos en InterIMAGE para hacer la segmentación. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 73 8.3. COMPOSICIÓN DE BANDAS Antes de utilizar InterIMAGE, debemos integrar en un mismo archivo las imágenes MDV Y NDVI. Con la herramienta Data Management \ Raster \ Raster Processing \ Composite Bands, como se muestra a continuación: Figura 109. Herramienta Composite Bands de ArcMap. Una vez realizada la operación, la imagen resultado tiene las siguientes características: Figura 110. Propiedades de la imagen formada por el MDV y el NDVI. Observamos que coinciden con el MDV y el NDVI. Exportamos la imagen en formato TIFF para poder trabajar con ella en InterIMAGE. Creamos un nuevo proyecto donde cargamos la imagen que integra tanto el MDV como el NDVI como se muestra a continuación. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 74 Figura 111. NDVI visualizado en InterIMAGE. También cargamos la ortofoto como elemento de apoyo. 8.4. REGLAS DE SEGMENTACIÓN Creamos las redes semánticas de forma individual, ya que el proceso que vamos a realizar se colapsa cuando introducimos más de una clase. Esto va influir de manera negativa, en cuanto a tiempo dedicado, pero como posteriormente se verá, merece la pena. Las redes semánticas tienen el siguiente formato: Figura 112. Redes Semánticas. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 75 El nodo Segmentation, realizará una segmentación con un parámetro de escala muy pequeño para poder identificar de una forma más precisa los segmentos que cumplan las condiciones establecidas. El valor de los parámetros será el siguiente: Figura 113. Parámetros del nodo Segmentation. Las condiciones establecidas (matorral, otros, suelo, edificios y árboles respectivamente) serán: Figura 114. Condiciones establecidas para cada clase. En este caso, utilizamos la opción Merge All, puesto que nos interesa unir los objetos para realizar una segmentación posterior. La segmentación posterior se realiza en los nodos Tree, Grass, Other, Soil y Buildings respectivamente, cuyos operadores son el TopDown TA_Baatz_Segmenter. Sus parámetros y sus características extraídas (promedio de alturas y perímetro) se observan en la página siguiente. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 76 Figura 115. Parámetros de los nodos hijos. Figura 116. Características exportadas . El resultado visualizado en ArcMap, con la ortofoto de fondo, se muestra en la figura siguiente: Figura 117. Rodales a partir de la integración del MDV y del NDVI. Delimitación de rodales a partir de la integración de datos LiDAR y de datos multiespectrales 77 Figura 118. Detalle de los rodales a partir de la integración del MDV y del NDVI. En esta última segmentación se puede apreciar todavía algún error que otro, como se observa en la imagen anterior hay zonas pertenecientes al suelo clasificadas como otros. Al igual que en las dos segmentaciones anteriores, vamos a realizar una evaluación mediante una nueva matriz de confusión, cuyos polígonos de muestra serán los mismos que en los apartados anteriores, teniendo en cuenta la clase otros como suelo. Tras realizar el recorte obtenemos el siguiente resultado: Figura 119. Polígonos resultado en la tercera segmentación. Determinación de parámetros estructurales de los rodales 84 En total genera 12.542 filas con los parámetros citados anteriormente. Los valores generados por la herramienta están multiplicados por un millón. Esto quiere decir que la elevación media en la primera fila de la tabla anterior (4055000), realmente significa que tiene una altura de 4.055 metros. Para poder continuar con el proceso, realizamos un Join con ArcMap entre la tabla resultado y la capa rodales utilizada como input en la operación de recorte. Unimos ambos datos mediante su campo común, es decir, su identificador. El resultado es el siguiente (las alturas ya están en metros): Figura 130. Tabla de atributos después de la unión visualizada en ArcMap. Con estos resultados en la tabla podemos observar detalles importantes a la hora de evaluar pero para que dicha evaluación pueda ser más visual y fácil de entender se van a crear planos individuales de cada variable y así podemos observar mejor las diferencias en la página siguiente. Determinación de parámetros estructurales de los rodales 85 1. ELEVACIÓN MEDIA. Figura 131. Plano elevación media. En el plano anterior se puede observar como las cubiertas relacionadas con el suelo desnudo y herbáceas, matorrales y árboles se representan con una escala de verdes (de más claro a más oscuro, respectivamente). Por otra parte las alturas superiores a 11 metros que corresponden a los edificios se representan con un rosa claro. En este caso al trabajar solo con elevaciones no podemos diferenciar entre los edificios y el resto de construcciones como en el caso de InterIMAGE. Se puede observar como tenemos los mismos problemas entre edificios y los árboles próximos a ellos debido a la geometría de los primeros (forma de escalón). Determinación de parámetros estructurales de los rodales 86 En cuanto a las cubiertas de vegetación se puede observar como la clasificación es muy similar a la de cubiertas realizada por InterIMAGE (debido a que los polígonos utilizados en InterIMAGE son los mismos que los utilizados en CloudMetrics) como se muestra a continuación: Figura 132. Comparativa entre la solución InterIMAGE (izqda.) y la de CloudMetrics (dcha.). En la imagen resultado de CloudMetrics (derecha) aparecen 4527 rodales de matorral y 4045 de árboles, mientras que, en la imagen resultado de InterIMAGE aparecen 4529 rodales de matorral y 4046 de árboles. Por lo tanto, podemos concluir que la evaluación de las cubiertas de vegetación ha sido correcta con una probabilidad, en cuanto a número de polígonos, mayor del 99 %. El plano resultado del cálculo de CloudMetrics cuyo parámetro es la elevación media se encuentra también en el apartado Anexos con su correcta maquetación. Determinación de parámetros estructurales de los rodales 87 2. PUNTOS POR RODAL. Figura 133. Plano de puntos por rodal. En el plano anterior se representa mediante una escala de rojos el número de puntos por rodal ordenados de mayor a menor (claro-oscuro). Se puede observar como los rodales con mayor superficie, es decir, los correspondientes a la clase suelo, son los que mayor número de puntos contienen. Esto es debido al siguiente razonamiento: si nuestra zona de estudio tiene una superficie de 184.730 m2 (406 columnas y 455 filas, cuyo píxel mide 1 m de lado, nos dan el resultado anterior) y el número total de puntos es de 112.500, podemos calcular el número de puntos por metro cuadrado mediante el cociente entre los puntos y la superficie, respectivamente. Determinación de parámetros estructurales de los rodales 88 Para poder calcular el número de puntos total en nuestra zona de estudio se podría actuar de varias maneras. Una de las formas sería mediante la herramienta CloudMetrics que nos devuelve el número de puntos por rodal. Simplemente sumando todos los valores de dicha columna nos devolverá el resultado anterior. Otra forma más compleja pero más profesional, sería mediante el uso del ejecutable de FUSION. Mediante los comandos Tools / Miscellaneous utilities / Examine LAS file headers el programa nos devuelve un breve resumen sobre la cabecera de la nube de puntos de nuestra zona como se muestra a continuación: Figura 134. Cabecera de la nube de puntos. Con los datos anteriores, llegamos a la conclusión que nuestra zona tiene 0.68 puntos/m2. Ese es el motivo por el cual los rodales cuya superficie es mayor (suelo) tienen mayor número de puntos mientras que los rodales de menor superficie (edificios) contienen menos puntos. Anteriormente se han remuestrado los puntos LIDAR a 1 metro puesto que este tamaño de píxel se adapta mejor a las condiciones del terreno. El plano resultado del cálculo de CloudMetrics cuyo parámetro es el número de puntos se encuentra también en el apartado Anexos con su correcta maquetación. Determinación de parámetros estructurales de los rodales 89 3. ELEVACIÓN PERCENTIL 75. Figura 135. Plano de elevación con percentil 75. El plano anterior, consiste en una representación de la elevación por rodal teniendo en cuenta el valor del 75 % de los puntos y no de la media como en el primer caso. Se puede observar como las cubiertas relacionadas con el suelo desnudo y herbáceas, matorrales y árboles se representan con una escala de verdes (de más claro a más oscuro, respectivamente). Por otra parte las alturas superiores a 11 metros que corresponden a los edificios se representan con un rosa claro. En este caso al trabajar solo con elevaciones no podemos diferenciar entre los edificios y el resto de construcciones como en el caso de InterIMAGE. Determinación de parámetros estructurales de los rodales 90 Se puede observar también, como en casos anteriores, que tenemos los mismos problemas entre edificios y los árboles próximos a ellos debido a la geometría de los primeros (forma de escalón y tamaño de los propios árboles). En cuanto a las cubiertas de vegetación vemos como la clasificación no es tan similar a la de cubiertas realizada por InterIMAGE como se muestra a continuación: Figura 136. Comparación entre el plano elevación y el plano elevación con percentil 75. Es curioso cómo se aprecian algunos rodales considerados, anteriormente, como matorrales se convierten en clase árboles. Esto es debido a que en esos rodales el 75 % de los puntos o más tiene una altura clasificada como árboles. En la imagen con percentil 75 (derecha) aparecen 4450 rodales de matorral y 4195 de árboles, mientras que, en la imagen resultado de InterIMAGE aparecían 4429 rodales de matorral y 4146 de árboles. Por lo tanto, podemos concluir que la evaluación de las cubiertas de vegetación ha sido correcta con una probabilidad, en cuanto a número de polígonos, mayor del 97 %. El plano resultado del cálculo de CloudMetrics cuyo parámetro es la elevación con percentil 75 se encuentra también en el apartado Anexos con su correcta maquetación. Determinación de parámetros estructurales de los rodales 91 4. DESVIACIÓN TÍPICA. Figura 137. Plano de desviaciones típicas. El plano anterior, consiste en una representación de la desviación típica por rodal utilizando para ello una escala de colores como se observa en la leyenda. Este plano es interesante para saber que rodales sufren cambios bruscos de altura. Por ejemplo, rodales considerados como edificios que contienen una minúscula superficie de suelo o rodales considerados como árboles que contienen una minúscula superficie de suelo y/o, en ocasiones, de matorral tendrán una desviación típica mayor (zonas representadas de color rojo) que el resto. Determinación de parámetros estructurales de los rodales 92 5. INTENSIDADES. Figura 138. Plano de intensidades. Como hemos dicho al principio del apartado, la herramienta CloudMetrics también opera con información relacionada con la intensidad. El plano anterior, consiste en una representación de la intensidad en Hz por rodal. La intensidad es una medida, recogida para cada punto, de la fuerza de retorno del pulso láser que genera el punto. Se basa en la reflectividad del objeto alcanzado por el pulso láser. Hay que tener en cuenta que la reflectividad es una función de la longitud de onda utilizada, que suele estar en el infrarrojo cercano. Si nos fijamos en la leyenda, se observa que el orden de la intensidad de las clases, de mayor a menor, es edificios, suelo, matorral y árboles Determinación de parámetros estructurales de los rodales 93 Figura 139. Plano detalle de las intensidades. Como representa el plano detalle anterior, podemos ver como el suelo (representado en un tono anaranjado) incluye zonas de tonos verdes que prueban la existencia de pequeñas herbáceas en la zona. La intensidad sirve de ayuda en la detección y extracción de entidades, en la clasificación de puntos LIDAR y como sustituta de imágenes aéreas cuando no hay ninguna disponible. Si los datos LIDAR incluyen valores de intensidad (como es el caso), se pueden crear imágenes a partir de ellos que parecen fotografías aéreas en blanco y negro o también se les puede asignar colores. 100 730700 730700 730800 730800 730900 730900 731000 731000 4360200 4360200 4360300 4360300 4360400 4360400 4360500 4360500 4360600 4360600 4360700 4360700 CLASIFICACIÓN POR CUBIERTAS (INTERIMAGE) rodales clases Suelo Otros Matorral Árboles Edificios AUTOR: RAFAEL LLORENS COMPANY ESCALA: 1:3.000 FECHA : 02/06/2016 SISTEMA DE COORDENADAS: ETRS89 HUSO 30 PROYECCIÓN UTM ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA Y TOPOGRÁFICA UNIVERSIDAD POLITÉCNICA DE VALENCIA El Pla del Garrofer Camí Vell de la Devesa Antic Tallafoc de la calle Platja de la Garrofera 730700 730700 730800 730800 730900 730900 731000 731000 4360200 4360200 4360300 4360300 4360400 4360400 4360500 4360500 4360600 4360600 4360700 4360700 CLASIFICACIÓN POR CUBIERTAS rodales clases Suelo Matorral Arboles Construcciones AUTOR: RAFAEL LLORENS COMPANY ESCALA: 1:3.000 FECHA : 02/06/2016 SISTEMA DE COORDENADAS: ETRS89 HUSO 30 PROYECCIÓN UTM ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA Y TOPOGRÁFICA UNIVERSIDAD POLITÉCNICA DE VALENCIA El Pla del Garrofer Camí Vell de la Devesa Antic Tallafoc de la calle Platja de la Garrofera 730700 730700 730800 730800 730900 730900 731000 731000 4360200 4360200 4360300 4360300 4360400 4360400 4360500 4360500 4360600 4360600 4360700 4360700 CLASIFICACIÓN POR ELEVACIÓN MEDIA rodales ele_media 0 - 0,2 m 0,2 - 3 m 3 - 11 m 11 - 29 m AUTOR: RAFAEL LLORENS COMPANY ESCALA: 1:3.000 FECHA : 02/06/2016 SISTEMA DE COORDENADAS: ETRS89 HUSO 30 PROYECCIÓN UTM ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA Y TOPOGRÁFICA UNIVERSIDAD POLITÉCNICA DE VALENCIA El Pla del Garrofer Camí Vell de la Devesa Antic Tallafoc de la calle Platja de la Garrofera 730700 730700 730800 730800 730900 730900 731000 731000 4360200 4360200 4360300 4360300 4360400 4360400 4360500 4360500 4360600 4360600 4360700 4360700 CLASIFICACIÓN POR ELEVACIÓN CON PERCENTIL 75 rodales elev_75 0 - 0,2 m 0,2 - 3 m 3 - 11 m 11 - 29 m AUTOR: RAFAEL LLORENS COMPANY ESCALA: 1:3.000 FECHA : 02/06/2016 SISTEMA DE COORDENADAS: ETRS89 HUSO 30 PROYECCIÓN UTM ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA Y TOPOGRÁFICA UNIVERSIDAD POLITÉCNICA DE VALENCIA El Pla del Garrofer Camí Vell de la Devesa Antic Tallafoc de la calle Platja de la Garrofera 730700 730700 730800 730800 730900 730900 731000 731000 4360200 4360200 4360300 4360300 4360400 4360400 4360500 4360500 4360600 4360600 4360700 4360700 CLASIFICACIÓN POR PUNTOS POR RODAL rodales Puntos 1 - 10 11 - 35 36 - 90 91 - 600 AUTOR: RAFAEL LLORENS COMPANY ESCALA: 1:3.000 FECHA : 02/06/2016 SISTEMA DE COORDENADAS: ETRS89 HUSO 30 PROYECCIÓN UTM ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA Y TOPOGRÁFICA UNIVERSIDAD POLITÉCNICA DE VALENCIA El Pla del Garrofer Camí Vell de la Devesa Antic Tallafoc de la calle Platja de la Garrofera