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) + gt y = (rows * gt) + gt X = np.arange(gt, x, gt) Y = np.arange(gt, y, gt) # 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)-1, axis = 1) X = np.delete(X, len(X)-1, axis = 1) Y = np.delete(Y, len(Y)-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()
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