python - Display a georeferenced DEM surface in 3D matplotlib

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

Answer

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

• the original DEM

• 
  

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

• 
  

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

!

The coloured grid (with matplotlib):

• 
  

  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)
  

With Visvis, you can make animations:

You can even plot the 3D contours:

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