Matplotlib 3D – Display a Georeferenced DEM Surface in 3D with Matplotlib

gdalmatplotlibnumpyosgeopython

I want to use a DEM file to generate a simulated terrain surface using matplotlib. But I do not know how to georeference the raster coordinates to a given CRS. Nor do I know how to express the georeferenced raster in a format suitable for use in a 3D matplotlib plot, for example as a numpy array.

Here is my python code so far:

import osgeo.gdal
dataset = osgeo.gdal.Open("MergedDEM")
gt = dataset.GetGeoTransform()

Best Answer

Matplotlib knows nothing about georeferenced surfaces, it only knows x,y,z coordinates. You can also use Visvis or Mayavi.

  • the original DEM

enter image description here

  • you must first extract the x,y, z coordinates of the grid ( a raster is a grid of pixels and with a DEM, the value of the pixel is the elevation, z) with osgeo.gdal. No script here because it is possible to find the solutions on Gis StackExchange or on the web.

  • after, you can plot the points in 3D

    from mpl_toolkits.mplot3d.axes3d import *
import matplotlib.pyplot as plt
from matplotlib import cm
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)  
plt.show()

enter image description here

  • and you must reconstruct a 3D grid (surface) with the function griddata of matplotlib (Delaunay)

    import numpy as np
from matplotlib.mlab import griddata
# craation of a 2D grid
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = griddata(x, y, z, xi, yi)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.show()

The grid (with visvis):

enter image description here!

The coloured grid (with matplotlib):

enter image description here

  • but you can also use others interpolation algorithms (Scipy thin-plate spline here) and drape colours

    import scipy as sp
import scipy.interpolate
# 2D grid construction
spline = sp.interpolate.Rbf(x,y,z,function='thin-plate')
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = spline(X,Y)

enter image description here

With Visvis, you can make animations:

enter image description here

You can even plot the 3D contours:

enter image description here

  • Comparison between a projected DEM (with GRASS GIS, x 1) and the equivalent non projected DEM (x,y,z, with Visvis x 1)

enter image description here enter image description here

Related Question