Tuesday, February 25, 2014

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.

14 comments:

  1. I really like that. Thanks for sharing!

    ReplyDelete
  2. thank you, but this is not working for me

    ReplyDelete
  3. here is my code(not working):
    bmp1 = vtk.vtkPNGReader()
    bmp1.SetFileName("157321.png")
    my_texture=vtk.vtkTexture()
    my_texture.interpolate=0
    my_texture.SetInput(0,bmp1.GetOutput())
    m=mlab.mesh(arrayx, arrayy, data.transpose(), color=(1,1,1))
    m.actor.enable_texture = True
    m.actor.tcoord_generator_mode = 'plane'
    m.actor.actor.texture = my_texture

    ReplyDelete
  4. can you help me please ?

    ReplyDelete
  5. Getting a memory error, any clues as to a cause for such?

    ReplyDelete
  6. I was able to get this working using Mayavi 4.4.3-9 by applying the patch suggested at the end of the following issue:

    https://github.com/enthought/mayavi/issues/211

    which is modifying one line in mayavi/components/actor.py. I also needed to change the way the Texture was created, using

    my_texture=tvtk.Texture(input_connection=bmp1.output_port, interpolate=0)

    instead of lines 17-19 in the original draping.py script.

    ReplyDelete
  7. This is great example! May I ask which mayavi and tvtk version you were using? I downloaded the exactly same data and script from Github, generated the 3D map successfully, but did not get the exact same result. The land jpg image did not expand to the whole surface, but only cover the center part. I saw similar problems and no one answered on stackoverflow, so I guess it may be related to the python library version. I will appreciate it if you can share your version!! Thank you!!

    ReplyDelete
    Replies
    1. I am having the same problem, but i can't find the solution anywhere.
      By chance did you find a solution to that issue?

      Delete
    2. This comment has been removed by the author.

      Delete
    3. Try Resizing your image to the DEM size using
      resize(DEM_height, DEM_width)
      That should work

      Delete
  8. Its very informative. You can now get hands on online services for 3D Architectural Renderings that too with 20% discount on all services. Visit here- | 3d architectural walkthrough services

    ReplyDelete
  9. I am trying to wrap an image (already read as an array and then converted with the code below) over a digital elevation model with the texture generated from the following piece of code,

    def texture_from_np(image):
    """
    Creates a tvtk.Texture object from a numpy image array
    """
    scalars = tvtk.UnsignedCharArray(name='RGB',
    number_of_components=3)
    # vtk uses column-major order (Fortan order), whereas
    # numpy array is in row-major order (C order) by default
    scalars.from_array(image.reshape(-1, 3, order='F'))
    tvtk_img = tvtk.ImageData()
    tvtk_img.dimensions = [image.shape[0], image.shape[1], 1]
    tvtk_img.point_data.scalars = scalars
    texture = tvtk.Texture()
    texture.set_input_data(tvtk_img)
    return texture

    However, the DEM appears to be blank. I tried to see if this is a contrast stretching issue by generating a random matrix of zeros and ones and wrap it over the DEM but it appears to be not as the scene is still blank.

    ReplyDelete
  10. Don't be overwhelmed, though, because you really don't have to know the exact dates of when these photos were taken. Just sort them by year, or even by decade if you cannot identify them all by year. 3dm8 website

    ReplyDelete