1

I want to add a picture over a surface that I made with mayavi and gdal.

I found a website (http://geoexamples.blogspot.ch/2014/02/3d-terrain-visualization-with-python.html) which explains how to do it but it doesn't work for me. I also saw that there were many unanswered questions about it on the forum.

Here is my program: Why isn't it working ?

import gdal
from mayavi import mlab
from tvtk.api import tvtk
import PIL
import numpy as np

# maido is the name of a mountain
# tipe is the name of a french school project

# 1) opening maido geotiff as an array
maido = gdal.Open('dem_maido_tipe.tif')
dem_maido = maido.ReadAsArray()

# 2) transformation of coordinates
columns = maido.RasterXSize
rows = maido.RasterYSize
gt = maido.GetGeoTransform()
ndv = maido.GetRasterBand(1).GetNoDataValue()

x = (columns * gt[1]) + gt[0]
y = (rows * gt[5]) + gt[3]

X = np.arange(gt[0], x, gt[1])
Y = np.arange(gt[3], y, gt[5])

# 3) creation of a simple grid without interpolation
X, Y = np.meshgrid(X, Y)

#Mayavi requires col, row ordering. GDAL reads in row, col (i.e y, x) order
dem_maido = np.rollaxis(dem_maido,0,2)
X = np.rollaxis(X,0,2)
Y = np.rollaxis(Y,0,2)

# 4) deleting the "no data" values
dem_maido = dem_maido.astype(np.float32)
dem_maido[dem_maido == ndv] = np.nan #if it's NaN, mayavi will interpolate

# delete the last column
dem_maido = np.delete(dem_maido, len(dem_maido)-1, axis = 0)
X = np.delete(X, len(X)-1, axis = 0)
Y = np.delete(Y, len(Y)-1, axis = 0)

# delete the last row
dem_maido = np.delete(dem_maido, len(dem_maido[0])-1, axis = 1)
X = np.delete(X, len(X[0])-1, axis = 1)
Y = np.delete(Y, len(Y[0])-1, axis = 1)

# 5) ploting the raster and draping a picture of maido over the surface
img_filename = "maido.jpg"
img = PIL.Image.open(img_filename)
bmp = tvtk.JPEGReader(file_name=img_filename)

t = tvtk.Texture()
t.interpolate=0
t.set_input(0, bmp.get_output())

mlab.figure(size=(1000, 1000))

surf = mlab.surf(X, Y, dem_maido, color=(1,1,1), warp_scale="auto")

surf.actor.actor.mapper.scalar_visibility = False
surf.actor.enable_texture = True
surf.actor.tcoord_generator_mode = 'plane'

mlab.show()

enter image description here

On the left you have what I get and on the right you have the picture I want to use.

Here is the GeoTiff file and the picture that I used (it is uploaded on my own google account don't worry): https://drive.google.com/drive/folders/0B7P95aWmH4DUbzJCaVZQM3Y0Y1U?usp=sharing

Your Answer

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

Browse other questions tagged or ask your own question.