3D terrain visualization with python and Mayavi2

I have always wanted to draw these 3D terrains like those in www.shadedrelief.com, which are amazing. But the examples were all using software I don't use, so I tried to do it with python.

The final result

As usual, you can get all the source code and data at my GitHub page.

 Getting the data

After trying different locations, I decided to use the mountain of Montserrat, close to Barcelona, since it has nice stone towers that are a good test for the DEM and the 3D visualization. An actual picture of the zone used is this one:

Montserrat monastery

The building is a good reference, since the stone only areas make the result testing much more difficult.
All the data has been downloaded from the ICGC servers:

• The DEM data was downloaded from the Vissir3 service, going to the section catàleg i descàrregues -> MDE 5x5. The file is named met5v10as0f0392Amr1r020.txt, but I cut a small part of it, to make mayavi2 work smoother using:
  
  gdalwarp -te 401620 4604246 403462 4605867 -s_srs EPSG:25831 -t_srs EPSG:25831 met5v10as0f0392Amr1r020.txt dem.tiff
• The picture to drap over the dem file can be downloaded using the WMS service given by the ICGC:
  
  http://geoserveis.icc.cat/icc_mapesbase/wms/service?REQUEST=GetMap&VERSION=1.1.0&SERVICE=WMS&SRS=EPSG:25831&BBOX=401620.0,4604246.0,403462.0,4605867.0&WIDTH=1403&HEIGHT=1146&LAYERS=orto5m&STYLES=&FORMAT=JPEG&BGCOLOR=0xFFFFFF&TRANSPARENT=TRUE&EXCEPTION=INIMAGE

 It's not as automatic as I would like, but if it's possible to download a DEM and the corresponding image, it's possible to create the 3D image.

Creating the image

First, let's plot the DEM file in 3D using mayavi2:

"""
Plotting the terrain DEM with Mayavi2
"""

from osgeo import gdal
from mayavi import mlab

ds = gdal.Open('dem.tiff')
data = ds.ReadAsArray()

mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

mlab.surf(data, warp_scale=0.2) 

mlab.show()

• In first place, we import as usual, gdal and numpy. Also the mlab library from mayavi, which lets set the mayavi canvas.
• The data is read, as usual, with the gdal ReadAsArray method. 
• The figure is created. This works like creating the Image object in the PIL library, creating the canvas where the data wil be drawn. In this case, the size is 640 x 800 pixels, making the figure bigger can affect the performance in some old computers. bgcolor sets the blue color as the background.
• The surf method will plot the surface. The input has to be a 2D numpy array, which is what we have.  
• The warp_scale argument sets the vertical scale. In this case, letting the default value (1?) creates a really exaggerated effect so its better to play a little to get a more realistic effect.
• The colors depend of the Z value at each point, and can be changed using the color or colormap option.
• The show() method makes the image to stay visible when running the example from a script. If you use ipython, you don't need this step.
• If you want to save the figure as a png, you can either use the icon in the mayavi window or call the method mlab.savefig('image_name')
• If you want to move the camera (change the prespective), you can use the roll/yaw/pitch methods:

  f = mlab.gcf()
  camera = f.scene.camera
  camera.yaw(45)

The plotted DEM

Now, let's put an aerial image over the 3D visualization:

"""
Draping an image over a terrain surface
"""
from osgeo import gdal
from tvtk.api import tvtk
from mayavi import mlab
import Image

ds = gdal.Open('dem.tiff')
data = ds.ReadAsArray()
im1 = Image.open("ortofoto.jpg")
im2 = im1.rotate(90)
im2.save("/tmp/ortofoto90.jpg")
bmp1 = tvtk.JPEGReader()
bmp1.file_name="/tmp/ortofoto90.jpg" #any jpeg file

my_texture=tvtk.Texture()
my_texture.interpolate=0
my_texture.set_input(0,bmp1.get_output())


mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

surf = mlab.surf(data, color=(1,1,1), warp_scale=0.2) 
surf.actor.enable_texture = True
surf.actor.tcoord_generator_mode = 'plane'
surf.actor.actor.texture = my_texture

mlab.show()

• The most important new import is tvtk. TVTK is a python api that allows to work with VTK objects. Actually, my knowledge of Mayavi2 is very limited, but I see TVTK as an extension.
• The DEM data is read the same way, using the ReadAsArray method.
• The aerial image, named ortofoto.jpg, is not in the correct orientation. It took me a lot of time to get what was happening. I rotate it opening the image with the PIL library and using the rotate method (lines 11 to 13)
• Then, the tvtk object with the texture is created, loading the image with a JPEGReader object, and assigning it to the Texture object (lines 14 to 19). 
• The figure and the 3D surface is created as in the other example (lines 22 and 24)
• Then, the surface is modified to show the image over it (lines 25 to 28). 

The final result

The result is correct, but the aerial image is a little moved from the place it should be so, since the terrain is really steep, part of the buildings are drawn on the stone walls! I edited the WMS coordinates a little so the result is slightly better. Anyway, the method is correct.

Links

• GitHub: The source code and example data
• Mayavi2: 3D scientific data visualization and plotting 
• This entry to a mailing list gave me the tip to create the example
• If you want to do more or less the same, but using JavaScript, Bjørn Sandvik posted this excellent example.
• ShadedRelief: Cool 3D and shaded relief examples.