I made a surface plot with matplotlib, mplot3d and gdal. Here is the code:
import gdal
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
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()
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)
# 4) deleting the "no data" values
# 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) plot the raster
fig, axes = plt.subplots(subplot_kw={'projection': '3d'})
surf = axes.plot_surface(X, Y, dem_maido, rstride=1, cstride=1, cmap=cm.gist_earth,linewidth=0, antialiased=False)
plt.colorbar(surf) # adding the colobar on the right
plt.show()
Here is the GeoTiff file that I used (it is uploaded on my own google account don't worry):
https://drive.google.com/file/d/0B7P95aWmH4DUQk9SbzhNUVNINGs/view
And here is what I get:
I am happy with the result but when I try to move it, it takes to much time.
What can I do to make it faster ?
Answer
mplot3d
is extremely slow because it uses software rendering. Watch your memory and CPU usage when you run that script, it will redline a CPU and use 1-2 GBs of RAM, just to render a tiny raster...
Downsampling will just reduce the quality/resolution of your plot. A better way to speed up your plot is to use an OpenGL capable 3D plot library, like Mayavi that will offload the rendering to your graphics card instead.
Mayavi can be difficult to install in Windows, the easiest way is to install a scientific python distribution like Anaconda.
Then you could use mayavi.mlab.surf
instead of axes.plot_surface
, see below for example. Rotating this interactively is instantaneous.
I leave it to you to add axes, change colour ramp and add colour bar (hint, you can play around manually with the plot, see 2nd screenshot below).
import gdal
#from mpl_toolkits.mplot3d import Axes3D
#from matplotlib import cm
#import matplotlib.pyplot as plt
from mayavi import mlab
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)
print (columns, rows, dem_maido.shape)
print (X.shape, Y.shape)
# 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) plot the raster
#fig, axes = plt.subplots(subplot_kw={'projection': '3d'})
#surf = axes.plot_surface(X, Y, dem_maido, rstride=1, cstride=1, cmap=cm.gist_earth,linewidth=0, antialiased=False)
#plt.colorbar(surf) # adding the colobar on the right
#plt.show()
surf = mlab.surf(X, Y, dem_maido, warp_scale="auto")
mlab.show()
And an example of what you can do:
No comments:
Post a Comment