Cómo extraer las coordenadas x, y, z de un ráster *.tif para visualización 3D con Python

En artículos anteriores (1, 2) se ha estado investigando acerca de la manera de cómo visualizar datos en 3D usando las librerías de representación gráfica de Python. En el artículo que me ha servido como referencia, se hace énfasis en la visualización en 3D de capas vectoriales (bien sean shapefiles de tipo 3D o shapefiles con atributo de altitud) o de rásteres en formato ASCII grid (*.asc). Sin embargo, se va a tratar de resumir todo el procedimiento con la consideración de un sólo tipo de formato, el ráster *.tif, porque, en principio, los vectoriales de punto pueden ser rasterizados con incorporación de la componente z como value y al hecho de que el formato *.asc, entre otras cosas, pesa el doble que el *.tif y no incorpora información acerca de la proyección.

La selección del formato *.tif se basa en el hecho de que descubrí la manera de cómo leer de manera eficiente el array que está embebido en ellos usando la librería GDAL de Python. El procedimiento se esboza en un artículo en el cual se determinan los valores promedios de un raster con diferentes formatos (*.jpg y *.tif) y tipo de datos. Por tanto, en lugar de retornar el promedio, la función respectiva se puede modificar ligeramente para que retorne en su lugar los tres arrays x, y, z. No obstante, en este caso particular (por premura), sólo voy a retornar unas posiciones de referencia de las componentes x, y (los valores proyectados pueden ser obtenidos fácilmente a través del procedimiento GetGeoTransform de GDAL). Esta parte del código va a ser acoplada antes de la visualización 3D mediante interpolación TIN 3D usando matplotlib. El código es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8

#Librerías
from osgeo import gdal
from os import system
import struct
import time

import numpy as np
from matplotlib.mlab import griddata
from mpl_toolkits.mplot3d.axes3d import *
from matplotlib import cm
import matplotlib.pyplot as plt

#Función que conforma los arreglos x,y,z
def getCoorXYZ(band):

	# fmttypes: Byte, UInt16, Int16, UInt32, Int32, Float32 y Float64
	fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

	print "filas = %d columnas = %d" % (band.YSize, band.XSize)

	BandType = gdal.GetDataTypeName(band.DataType)

	print "Tipo de datos = ", BandType

	totHeight = 0
	
	x = []
	y_ = []
	z = []

	inc_x = 0

	for y in range(band.YSize):
	 
		scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
		values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
	 		
		for value in values:
			z.append(value)
			inc_x += 1
			y_.append(inc_x)
			x.append(y+1)			
		
		inc_x = 0
	     
	return x, y_, z

#Comienzo del programa

system("clear")
	 
nameraster = str(raw_input("Nombre del ráster = ? "))

start = time.time()
	 
dataset = gdal.Open(nameraster)
band = dataset.GetRasterBand(1)

print "Procesando %s" % nameraster

x,y,z = getCoorXYZ(band)

# construcción de la grilla 2D
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
  
# interpolación
Z = griddata(x, y, z, xi, yi)
  
#Visualización con Matplotlib
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.plot
  
end = time.time()

time_tot = end - start

print "Tiempo total = %.4f s" % time_tot	 

plt.show() #Para que la ventana de visualización permanezca estática

El ráster que usé para la visualización fue el que obtuve aquí:

Cálculo de la longitud del terreno (λ de la RUSLE) según enfoque CALSITE mediante un script de Python con GRASS y GDAL

por lo que la componente z es, en este caso, la longitud del terreno (λ de la RUSLE). Esta es la visualización que se produce al ejecutar el script (rotado convenientemente para obtener una buena perspectiva):

extraccion

y la información complementaria que se obtiene, en cónsola de bash, es la siguiente:

Nombre del ráster = ? salida2.tif
Procesando salida2.tif
filas = 50 columnas = 50
Tipo de datos =  Float64
Tiempo total = 0.6635 s

ejecutándose en 663.5 milésimas de segundo.

Para saber si el código se ejecutó de manera adecuada voy a rotar el mapa 2D para ver si es equivalente al 3D.

extracción2

En este caso lo son. Sin embargo, el código que había producido con anterioridad a éste invertía las posiciones x, y por lo que producía la imagen especular del mapa verdadero. Está es “la prueba del ácido” para verificar que todo está funcionando, aparentemente, de la manera adecuada

Esta entrada fue publicada en Código Python, SIG, Software Libre. Guarda el enlace permanente.

4 respuestas a Cómo extraer las coordenadas x, y, z de un ráster *.tif para visualización 3D con Python

  1. Pingback: Visualización 3D con Python (gdal, numpy y matplotlib) en Windows | El Blog de José Guerrero

  2. Pingback: Cómo extraer las coordenadas x, y, z de ...

  3. Pingback: Paquete ‘rgl’ para la visualización 3D en lenguaje R | El Blog de José Guerrero

  4. Pingback: Subregiones ráster usando GDAL en ambiente de PyQGIS | El Blog de José Guerrero

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s