From Display a georeferenced DEM surface in 3D matplotlib
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)
