Hoy vamos a utilizar la librería matplotlib para obtener con Python un gráfico tridimensional de un modelo digital del terreno (MDT). Esta entrada se basa en el ejemplo que nos da matplotlib de gráficos 3D, pero previamente vamos a preparar un fichero .npz a partir de un ráster de entrada.
Obtener un MDT
Si disponemos de un ráster MDT podemos saltarnos este paso. En caso de que no sea así, nos dirigimos al centro de descargas del CNIG y en el apartado «Modelos Digitales de Elevaciones» seleccionamos la precisión y luego en el mapa el fichero de la zona que queramos. Esto nos descarga un fichero .asc que podemos cargar en QGIS para ver el aspecto que tiene.
Convertir nuestro MDT a formato .npz
Antes de nada cargaremos los siguientes módulos en Python:
from matplotlib import cbook from matplotlib import cm from matplotlib.colors import LightSource from mpl_toolkits.mplot3d import axes3d, Axes3D import matplotlib.pyplot as plt import numpy as np import gdal
El formato .npz (el que se usa de entrada en el ejemplo) es un fichero de numpy que contiene diferentes arrays. Vamos a guardar los valores que van a ser relevantes para la generación del gráfico y para ello necesitamos obtenerlos del ráster. Aquí es donde Python entra en juego para obtener la elevación y las esquinas coordenadas de las esquinas que delimitan los datos. Para ello usaremos el siguiente código:
#Cargamos la ruta del ráster en la variable "dsm" antes image = gdal.Open(dsm) band = image.GetRasterBand(1) gt = image.GetGeoTransform() xsize = gt[1] #Leemos los valores del MDT como un array de dimensiones x, y (x arrays de y elementos) z = band.ReadAsArray() nrows, ncols = z.shape xmin = gt[0] xmax = xmin + (ncols * xsize) ymin = gt[2] #Conocemos la resolución en "y", es igual a la resolución en "x" ymax = ymin + (nrows * xsize)
Ahora ya tenemos almacenados en las variables z, xmin, xmax, ymin e ymax los valores que nos interesan. Vamos a guardarlos en un fichero .npz mediante numpy y estableciendo parejas clave-valor donde elevation serán los valores que tenemos en la variable «z».
#Hemos importado numpy como np #npz es el nombre de la variable que contiene la ruta para el fichero npz np.savez(npz, elevation = z, xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
Ya tenemos un fichero .npz a partir de nuestro MDT. Y ahora procederemos a generar el gráfico.
Uso de matplotlib para generar el gráfico tridimensional
Mediante las siguientes líneas abriremos el archivo y le daremos forma a los datos.
with cbook.get_sample_data(npz) as file, \ np.load(file) as dem: z = dem['elevation'] nrows, ncols = z.shape x = np.linspace(dem['xmin'], dem['xmax'], ncols) y = np.linspace(dem['ymin'], dem['ymax'], nrows) x, y = np.meshgrid(x, y)
Es necesario ahora definir una región de datos que es la que se va a imprimir. Si los datos no son muy extensos se puede usar np.s_[0: , 0: ], aunque el método también sirve para recortar valores anómalos que podamos encontrar en los bordes de la imagen.
region = np.s_[10:-10, 10:-10] x, y, z = x[region], y[region], z[region]
Ahora procedemos a configurar el espacio:
fig = plt.figure() ax = Axes3D(fig) ls = LightSource(270, 45) rgb = ls.shade(z, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft') surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=rgb, linewidth=0, antialiased=False, shade=False) #Invierte la dirección del eje x ax.invert_xaxis() #Con plt.show() abrimos la ventana que mostrará el gráfico plt.show()
Hasta aquí el código nos devuelve la imagen que encabeza el artículo que representa el Puerto del León en Málaga. Pero podemos darle otra rampa de color mediante:
N = 256 vals = np.ones((N, 4)) #El primer número refleja el color en la escala de 256 valores de RGB para el canal 0=R, 1=G, 2=B vals[:, 0] = np.linspace(10/256, 1, N) vals[:, 1] = np.linspace(35/256, 1, N) vals[:, 2] = np.linspace(5/256, 1, N) newcmp = ListedColormap(vals)
Y obtenemos el siguiente resultado:
Como ves Python es una herramienta potente (y no solo en el ámbito de los SIG) si quieres aprender más apúntate a alguno de nuestros cursos y recibirás formación de calidad impartida por profesionales. Para obtener más información sobre los trámites puedes escribir a formacion@tycgis.com.
Buenas tardes, en que consola se realizó el ejercicio, gracias.
Buenos días John,
puedes realizar el ejercicio en la consola de Python de QGIS o ejecutarlo en PowerShell o cmd… El único requisito imprescindible es que tengas instaladas todas las librerías necesarias. Tanto en QGIS como en cualquiera de los otros matplotlib se encargará de abrir la ventanita con el resultado.
Un saludo.
hola, mira tengo una pregunta, podría mostrar capas, como por ejemplo una vía encima del modelo DEM? Gracias
Buenos días, puedes mostrar todos los datos que añadas al espacio de impresión en matplotlib. Sin embargo, ten en cuenta que según qué usos puede que esta metodología no sea la idónea y existen otras herramientas.
Buenas tardes estimado una consulta Ud. cree que con esto pueda hacer un inventario de fallas geologicas con imagenes 3d
Buenos días Cristian, puede representar cualquier tipo de datos que pueda convertir a un array npz. Un ráster debería servir. Solo debe tener en cuenta que si la vista 3D no es manipulable para el usuario final puede que no sea la mejor representación de la información ya que una imagen plana se interpreta de un solo vistazo mientras que una vista isométrica en ocasiones no es una representación completa de toda la información.