Display a georeferenced DEM surface in 3D matplotlib
Question: 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: You can use the normal `plot_surface` method from matplotlib. Because it needs
a X and Y array, its already plotted with the right coordinates. I always find
it hard to make nice looking 3D plots, so the visual aspects can certainly be
improved. :)

    
    
    import gdal
    from mpl_toolkits.mplot3d import Axes3D
    
    dem = gdal.Open('gmted_small.tif')
    gt  = dem.GetGeoTransform()
    dem = dem.ReadAsArray()
    
    fig, ax = plt.subplots(figsize=(16,8), subplot_kw={'projection': '3d'})
    
    xres = gt[1]
    yres = gt[5]
    
    X = np.arange(gt[0], gt[0] + dem.shape[1]*xres, xres)
    Y = np.arange(gt[3], gt[3] + dem.shape[0]*yres, yres)
    
    X, Y = np.meshgrid(X, Y)
    
    surf = ax.plot_surface(X,Y,dem, rstride=1, cstride=1, cmap=plt.cm.RdYlBu_r, vmin=0, vmax=4000, linewidth=0, antialiased=True)
    
    ax.set_zlim(0, 60000) # to make it stand out less
    ax.view_init(60,-105)
    
    fig.colorbar(surf, shrink=0.4, aspect=20)
    

![enter image description here](http://i.stack.imgur.com/UCktv.png)

