Visualización 3D con Python (gdal, numpy y matplotlib) en Windows

Hace unos días publiqué unos posts relacionados con la visualización 3D en GNU/Linux usando Python y las librerías gdal, numpy, matplotlib y visvis (1, 2, 3). Hoy voy a explorar la posibilidad de hacerlo en Windows. Para ello hay que instalar conjuntamente con Python las librerías gdal, numpy y matplotlib. A diferencia de GNU/Linux, la instalación puede producir “dolores de cabeza”, sobre todo en los poco experimentados, porque hay que tener cuidado en que las versiones sean compatible; algo que en GNU/Linux se gestiona automáticamente a través de repositorios (sobre todo para los que tienen una Debian o Debian derivada). Como tengo la versión 2.7 de Python, la de numpy fue esta:

numpy-1.5.1-win32-superpack-python2.7.exe

La de matplotlib la bajé de acá:

https://github.com/matplotlib/matplotlib/downloads

y para instalar gdal hay que seguir meticulosamente este tutorial:

http://cartometric.com/blog/2011/10/17/install-gdal-on-windows/

Una vez instaladas y probadas en el intérprete que se cargaban las librerías, ejecuté este código:

#Libreri­as
from osgeo import gdal
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

#Funcion 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

nameraster = str(raw_input("Nombre del raster = ? "))

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

print "Procesando %s" % nameraster

x,y,z = getCoorXYZ(band)

# construccion de la grilla 2D
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
 
# interpolacion TIN
Z = griddata(x, y, z, xi, yi) #requiere 5 argumentos
 
#Visualizacion 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)

#Etiquetas de los ejes X, Y
plt.xlabel('X', fontsize=14, color='blue')
plt.ylabel('Y', fontsize=14, color='blue')

plt.plot
  
end = time.time()

time_tot = end - start

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

plt.show() #Para que la ventana de visualizacion permanezca estatica

que produjo este resultado:

python3_D

El código, más el *.tif de muestra, se pueden bajar de aquí:

python_dem.zip

Esta entrada fue publicada en Código C++, SIG, Suelos y Aguas. Guarda el enlace permanente.

4 respuestas a Visualización 3D con Python (gdal, numpy y matplotlib) en Windows

  1. david arturomendez matute dijo:

    bueno gracias

  2. Pingback: Interpolación spline cubica con gdalwarp de mapas de precipitación obtenidos por interpolación TIN | El Blog de José Guerrero

  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