Superficie de interpolación TIN 3D con los módulos numpy, matplotlib y visvis de Python

En el artículo anterior se había considerado la manera de obtener una superficie de interpolación spline 3D con los módulos scipy y visvis de Python. En este nuevo artículo se considera ahora una superficie de interpolación TIN 3D pero haciendo principalmente énfasis en las prestaciones de matplotlib de Python (disponible en los repositorios de Debian).

En el post anterior prácticamente se comenzó por el final del artículo de Martin Laloux porque considero que los aspectos más críticos están relacionados con la visualización 3D. Más adelante trataré la recuperación de las coordenadas.

Por otra parte, cuando comencé a analizar la porción mínima de código que me llevaría a una posible superficie de interpolación TIN 3D me arrojó error cuando tomé los arrays x, y, z del post anterior. La función griddata de matplotlib señalaba la imposibilidad de ejecutar la triangulación de Delaunay para esa particular disposición de puntos “inventados”. Era obvio porque estos no representaban una superficie sino una curva (cuya proyección en el plano z = 0 es una recta) en el espacio tridimensional (que para la interpolación spline si funciona). Por tanto, lo único que tuve que hacer fue “extender” los puntos y con sólo 16 se produjo un resultado con el código que se presenta a continuación:

#!/usr/bin/env python
#-* coding: utf-8
 
#Librerías
import numpy as np
from matplotlib.mlab import griddata
from mpl_toolkits.mplot3d.axes3d import *
from matplotlib import cm
import matplotlib.pyplot as plt
 
#Arreglos de puntos: ahora son 16!
x=[0,1,2,3,4,5,0,0,0,0,0,1,2,3,4,5]
y=[0,0,0,0,0,0,1,2,3,4,5,1,2,3,4,5]
z=[1,5,1,3,2,1,1,7,1,4,2,1,2,2,2,1]
 
# 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
 
plt.show() #Para que la ventana de visualización permanezca estática

La visualización fue la siguiente:

spline2

Para obtener lo mismo, pero visualizado como una malla gris con visvis, se tiene este código (donde el valor de z ha sido atenuado a la mitad con daspect)

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

#Librerías
import numpy as np
from matplotlib.mlab import griddata
import visvis

#Arreglos de puntos: ahora son 16!
x=[0,1,2,3,4,5,0,0,0,0,0,1,2,3,4,5]
y=[0,0,0,0,0,0,1,2,3,4,5,1,2,3,4,5]
z=[1,5,1,3,2,1,1,7,1,4,2,1,2,2,2,1]

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

f = visvis.gca()
m = visvis.grid(xi,yi,Z) 
#m.colormap = visvis.CM_JET #Descomente para un mapa de color
f.axis.visible = False
f.daspect = 1,1,0.5 # z x 0.5; cambie el valor de z para atenuar o exagerar

# Necesario para que la vista no se cierre
app = visvis.use() 
app.Run()

El resultado de la ejecución del código anterior es el siguiente:

spline3

Finalmente, empleando visvis y el mismo factor de atenuación para z, se propone este código:

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

#Librerías
import numpy as np
from matplotlib.mlab import griddata
import visvis

#Arreglos de puntos: ahora son 16!
x=[0,1,2,3,4,5,0,0,0,0,0,1,2,3,4,5]
y=[0,0,0,0,0,0,1,2,3,4,5,1,2,3,4,5]
z=[1,5,1,3,2,1,1,7,1,4,2,1,2,2,2,1]

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

f = visvis.gca()
m = visvis.surf(xi,yi,Z) 
m.colormap = visvis.CM_JET 
f.axis.visible = False
f.daspect = 1,1,0.5 # z x 0.5; cambie el valor de z para atenuar o exagerar

# Necesario para que la vista no se cierre
app = visvis.use() 
app.Run()

que resulta en un relieve sombreado con una tabla de color CM_JET:

spline4

Ahora sólo falta la recuperación de coordenadas; tema de próximos posts.

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

3 respuestas a Superficie de interpolación TIN 3D con los módulos numpy, matplotlib y visvis de Python

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

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

  3. Pingback: Paquete ‘rgl’ para la visualización 3D en lenguaje R | 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