MDT en 3 dimensiones con Python

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.

(7 votos, promedio: 5,00 de 5)

Cargando...

Formación de calidad impartida por profesionales