[GIS] DEM plot with matplotlib is too slow

gdalgeotiff-tiffmatplotlibpython

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:

enter image description here

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 ?

Best Answer

Just to start, you could plot only a subset of the data. Maybe downsample by a factor of 4. Of course, you will need to think about resample method which would probably be bilinear in this case. What about something like the answer from this question here: https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array

def rebin(a, shape):
   sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
   return a.reshape(sh).mean(-1).mean(1)

# 1) opening maido geotiff as an array
maido = gdal.Open('dem_maido_tipe.tif')
dem_maido = maido.ReadAsArray()
resampled_dem = rebin(dem_maido, (y/4, x/4))
Related Question