Having trouble to get a nice plot of a DEM

Question

With the explanations that I got on this post (Plot a surface with a DEM and mplot3d) I was able to plot a DEM (that I made with QGIS) with matplotlib on Python.

However I am still having trouble to get a nice plot.

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) plot the raster
fig, axes = plt.subplots(subplot_kw={'projection': '3d'})
axes.plot_surface(X, Y, dem_maido, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)

Here is what I get:

Here is the GeoTiff file that I used (it is uploaded on my own google account don't worry): https://drive.google.com/open?id=0B7P95aWmH4DUQk9SbzhNUVNINGs

I am getting huge negative values and I don't know why.

EDIT: I deleted the negative values in the last row and the last column

print "\n------------------\n"
print dem_maido

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)

print "\n------------------\n"
print dem_maido

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)

print "\n------------------\n"
print dem_maido

gene · Accepted Answer · 2017-01-03T16:11:42.767

You must take into account the nodata values:

maido_bnd = maido.GetRasterBand(1)
print.maido_bnd = maido.band.GetNoDataValue()
-32768.0

If I compute the min and max values of the DEM

print dem_maido.min(),dem_maido.max()
-32768 2885

Therefore, the elevations are from -32768 to 2885 and the result is your figure with the values < 0 = nodata values.

And your last correction (deleting the negative values in the last row and the last column) is only valid in this particular case (with scatter3D)

If I try your script with

ax.set_zlim(0, 2900)

We can see the problem caused by these nodata value

In a more general case you need to extract only the valid x,y, z coordinates

from osgeo import gdal
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
maido = gdal.Open('dem_maido_tipe.tif')
maido_bnd = maido.GetRasterBand(1)
no_data=maido_bnd.GetNoDataValue()
dem_maido = maido_bnd.ReadAsArray()
gtr = maido.GetGeoTransform()
y, x = np.where(dem_maido != no_data)
z = np.extract(dem_maido != no_data, dem_maido)
gtr_x = np.array(gtr[0] + x * gtr[1] + (y* gtr[2])
gtr_y =  np.array(gtr[3] + x * gtr[4] + y * gtr[5])
print z.min(), z.max()
0 2885

I plot the x,y,z values of the new DEM (elevations from 0 to 2885)

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)  
plt.show()

The result is different and correct.

If you want a surface, you need to interpolate a new grid with the new values

from matplotlib.mlab import griddata
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
from matplotlib.mlab import griddata # Delaunay
Z = griddata(x, y, z, xi, yi)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=True)
plt.show()

Having trouble to get a nice plot of a DEM

1 Answers1