0

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:

enter image description here

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

enter image description here

2

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)

enter image description here

If I try your script with

ax.set_zlim(0, 2900)

We can see the problem caused by these nodata value

enter image description here

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.

enter image description here

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.