-2

I want to plot a surface on Python with a DEM that I made on QGIS.

I found 2 posts on this topic which were published by the same person:

  1. GIS Stackexchange: on this post they show us how to obtain a 3D surface from the data of elevations recovered from the DEM but they do not explain how they recovered the data;
  2. Stackoverflow: and on this post they show us how to recover the elevation data using the DEM and then they get a 3D surface

I understand everything explained on GIS Stackexchange because I am used to using matplotlib but I do not understand anything that is done in the Stackoverflow code (because I never used GDAL on Python).

Can you help me understand this code please ?

enter image description here

4
  • 1
    Its a bit difficult to say, since we don't know what you know. Can you ask a specific question? Have you looked at the documentation for the various methods being called?
    – BradHards
    Dec 20, 2016 at 9:12
  • I know how to plot a surface with the interpolation method when I have a list of X,Y and Z but I a bit ignorant for the rest (I hope my English is understandable by the way). I looked this documentation to understand the functions used: gdal.org/classGDALDataset.html. Dec 20, 2016 at 9:19
  • 1
    Can you try asking a specific question? Right now, its mostly "help me" which isn't that easy.
    – BradHards
    Dec 20, 2016 at 10:30
  • 1
    Please include any code as text rather than pictures.
    – PolyGeo
    Dec 20, 2016 at 10:37

1 Answer 1

4

1) In the first case, we have only the points and we need to construct a grid (regular or irregular) by interpolation.

2) In the second case, the raster is already a grid (regular) but Matplotlib only knows the x,y dimension of the image that you can plot

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from osgeo import gdal
test = gdal.Open('a.tif')
dem = test.ReadAsArray()
cols = test.RasterXSize
rows = test.RasterYSize
X = np.arange(0, cols)
Y = np.arange(0, rows)
X, Y = np.meshgrid(X, Y)
fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
surf = ax.plot_surface(X, Y, dem, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)

enter image description here

But matplotlib knows nothing about the georeferenced grid (look at Python affine transforms)

The typical geospatial coordinate reference system is defined on a cartesian plane with the 0,0 origin in the bottom left and X and Y increasing as you go up and to the right. But raster data, coming from its image processing origins, uses a different referencing system to access pixels. We refer to rows and columns with the 0,0 origin in the upper left and rows increase and you move down while the columns increase as you go right. Still a cartesian plane but not the same one.

The Gdal geotransform parameters allow to convert the coordinates between the two

from osgeo import gdal
test = gdal.Open('a.tif')
gt = test.GetGeoTransform()
print gt
(204073.37, 10.0, 0.0, 89817.66, 0.0, -10.0)

0 = x-coordinate of the upper-left corner of the upper-left pixel
1 = width of a pixel
2 = row rotation (typically zero)
3 = y-coordinate of the of the upper-left corner of the upper-left pixel
4 = column rotation (typically zero)
5 = height of a pixel (typically negative)

To transform the coordinates

# col, row to x, y
x = (col * gt[1]) + gt[0]
y = (row * gt[5]) + gt[3]

Therefore in the second solution

 from mpl_toolkits.mplot3d import Axes3D
 import matplotlib.pyplot as plt
 from matplotlib import cm
 import numpy as np
 dem = test.ReadAsArray()
 # transformation of coordinates
 xres = gt[1]
 yres = gt[5]
 X = np.arange(gt[0], gt[0] + dem.shape[1]*xres, xres)
 Y = np.arange(gt[3], gt[3] + dem.shape[0]*yres, yres)
 # creation of a simple grid without interpolation
 X, Y = np.meshgrid(X, Y)
 # plot the raster
 surf = ax.plot_surface(X, Y, dem, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)

enter image description here

You can also use Mayavi to visualize the 3D surface

enter image description here

And visualize the intersection of two surfaces (grids)

enter image description here

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