6

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

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

  • I am new in matplotlib, I converted my raster file to csv file in x y z format. How can I use my csv file in first your code? – nickan Oct 31 '17 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.