3

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()
Saullo G. P. Castro
  • 56,802
  • 26
  • 179
  • 234
EricVonB
  • 238
  • 3
  • 11
  • why do you need to georeference your DEM? – ala Jul 18 '13 at 09:35
  • I need to produce a surface in 3D matplotlib such that the Z component is the height as given in the DEM and the X and the Y components are the eastings and the northings from the DEM. – EricVonB Jul 18 '13 at 18:26

1 Answers1

4

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

Richard
  • 56,349
  • 34
  • 180
  • 251
Rutger Kassies
  • 61,630
  • 17
  • 112
  • 97
  • The code doesn't seem to run as it is posted. The imports for numpy and plt are missing and array dims don't match – Stein Dec 30 '18 at 17:50