# Having trouble to get a nice plot of a DEM

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')

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

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)

# 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)-1, axis = 1)
X = np.delete(X, len(X)-1, axis = 1)
Y = np.delete(Y, len(Y)-1, axis = 1)

print "\n------------------\n"
print dem_maido
`````` 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()
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 + x * gtr + (y* gtr)
gtr_y =  np.array(gtr + x * gtr + y * gtr)
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()
`````` 