23 noviembre 2012

Modelo Digital del Terreno (MDT) partiendo de datos LiDAR del PNOA

En este post vamos a descbir como generar un MDT en formato raster, de 5 m de tamaño de pixel, partiendo de de los ficheros disponibles en la web del CNIG procedentes de datos LiDAR del Plan Nacional de Ortofotografía Aérea (PNOA).

Para descargarnos los ficheros correspondientes a nuestra zona de estudio, tenemos que entrar en la web del centro de descargas del Centro Nacional de Información Geográfica(CNIG).

Dependiendo del uso que vayamos a hacer de la información (comercial o no comercial) tendremos que acceder por diferentes enlaces. En este caso como no se trata de ningún uso comercial, para descargar de la información debemos registrarnos como usuario para poder acceder a la zona de descarga de datos.

Bien por búsqueda avanzada o en el visor, accedemos a la descarga del Producto deseado. Buscamos el producto que nos interese, en este caso el Modelo Digital del Terreno - MDT05/MDT05-LIDAR, y en División Administrativa, seleccionamos por ejemplo la Provincia de Las Palmas.

Revisando la cuadrícula del Mapa Topográfico Nacional, escala 1:50.000 - MTN50, vemos que la Isla de gran Canaria se encuentra comprendida en las hojas número 1098, 1003, 1104, 1106, 1107. Pero nos falta la parte sur de la isla.

Si nos vamos al visor, vemos que aparece la hoja que nos falta en la parte sur de la isla, la 1109. Misterios de la geografía, debe de ser que San Borondón está cerca y por eso aparecen y desaparecen las hojas cartográficas.

Desde el centro de descargas buscamos las hojas que nos interesen y procedemos a descargalas, tras registrarnos como usuario y aceptar el acuerdo de licencia.

El formato en que descargaremos los archivos es ASCII Grid (.asc) . Este tipo de formato no es más que un fichero de texto que contiene una cabecera, donde se indican la número de columnas y de filas, coordenadas del pixel central, tamaño de pixel y valores de pixel sin datos (NoData)y una cuadrícula de valores, en este caso cada uno de los valores se corresponde con la elevación sobre el nivel del mar (Z ortométrica) de cada uno de los pixeles.

Para este ejemplo vamos a utilizar el software comercial ArcGIS de ESRI, en su versión de prueba de 60 días de duración, que se puede descargar haciendo clic aquí , tras crearse una Cuenta Global de Esri.

Una vez arrancamos ArcMap, vamos a añadir capas y navegamos hasta la carpeta donde hayamos descargados los ficheros .asc. Directamente que ArcMap ya los reconoce como una capa raster que podemos añadir a nuestro mapa.

Pero la gestión que hace de esta información es muy lenta, por lo que vamos a convertir a formato raster cada uno de los archivos, y crear un mosaico, generando una única capa raster con el MDT para la isla de Gran Canaria.

El primer paso es convertir de ASCII Grid (.asc) a formato raster (.tif). Para ello dentro de ArcToolbox, utilizaremos la herramienta Conversion Tools ♣ To Raster ♣ ASCII to Raster

Como son varios archivos y evitar repetir el proceso de conversión para cada uno de ellos, abriremos la herramienta en modo múltiple, para ello nos situamos sobre ella y hacemos clic con el botón derecho y seleccionamos Batch.

Una vez abierta la herramienta en modo múltiple, rellenamos cada uno de los campos de la tabla:
  • Input: cada uno de los ficheros .asc
  • Output raster: dejamos la Geodatabase por defecto (Default.gdb) y le damos nombre al fichero raster de salida, sin extensión. El nombre del archivo de salida no puede comenzar por caracteres numéricos
  • Output data: elegimos Float (queremos formato decimal para los valores de elevación) decimales)
Una vez completado el primer registro con el botón Add row (+), añadimos nuevos registros, tantos como ficheros deseamos procesar, en este caso cinco más. Una vez añadidos, cada uno de los registros del primera fila (Input / Output raster / Output data), hacemos clic con el botón derecho del ratón y seleccionamos, Fill, de esta manera los campo se rellenan automáticamente de forma análoga al primero. Ya solo tenemos que cambiar la parte final de las rutas; la parte correspondiente a los nombres ficheros de entrada y de los de salida. Una vez completo pulsamos aceptar y el resultado se añadirá directamente a nuestro mapa. Si esto no ocurre, desde el Menú Geoprocessing, Geoprocessing Options, en la parte inferior del cuadro de diálogo elegimos Add results of geoprocessing operations to the display, para seleccionar que los resultados de las operaciones de geoprocesamiento se carguen automaticamente en el proyecto.


En la imagen anterior, podemos observar el resultado de la conversión de los 6 archivos .asc a formato raster.

Utilizando la barra de herramientas Effects, y la herramienta Swipe, podemos comprobar cómo cada una de las imágenes contiguas tienen una banda de solape, que evita problemas de efecto borde o frontera entre escenas.

El último paso es crear un único raster que contenga a todas las escenas, para ello vamos a ArcToolbox, Data Management Tools ♣ Raster ♣ Raster Dataset ♣ Mosaic to New raster

Si descargamos los recursos, que nos dan información sobre los archivos descargados podemos leer en el fichero MDT05_LEEME.txt, entre otra la siguiente información:
  • MDT05: Modelo Digital del Terreno con paso de malla 5m
  • Sistema de Referencia Geodésico ETRS89 (REGCAN95 en Canarias, compatible con ETRS89/WGS84)
  • Proyección Cartográfica UTM en el huso correspondiente
Por lo tanto ya sabemos el tamaño de pixel de salida (5 m), el sistema de referencia geodésico REGCAN95, con la proyección UTM y la zona 28N.
Seleccionamos el tipo de pixel de salida, 32_BIT_FLOAT , el número de bandas, que en este caso es una y definimos de donde se obtendrán los datos en las zonas de solape, dejando la opción por defecto LAST.

Eligiendo una paleta adecuada (de un tono claro para valores y que se oscurezca a medida que aumenta el valor de altitud), coloreando entre los valores mínimos y máximo, eliminado los valores de fondo y dando efecto de sombreado este es el resultado final para el MDT de la isla de Gran Canaria en formato .tif derivado de datos LiDAR del PNOA y con un tamaño de pixel de 5m.

Sin más que un rápido análisis visual podemos detectar en la zona sureste de la isla alguna zona sin datos. Estas zonas sin datos provienen de los datos LiDAR de partida y no con la metodología seguida. Puede ser que la zona en cuestión estuviese cubierta de agua durante la adquisición de los datos. El LiDAR en zonas de agua es absorbido en su totalidad y no produce ningún retorno y se producen esas zonas sin datos (NoData). Si nos interesa rellenarlas podríamos aplicar algún tipo de filtro sobre esa zona.

1 comentario:

Unknown dijo...
Este comentario ha sido eliminado por un administrador del blog.

Publicar un comentario