Отобразить поверхность матрицы высот с географической привязкой в ​​3D Matplotlib

Question

Я хочу использовать файл DEM для генерации имитируемой поверхности рельефа с помощью matplotlib. Но я не знаю, как привязать растровые координаты к заданному CRS. Также я не знаю, как выразить растр с географической привязкой в ​​формате, подходящем для использования на трехмерном графике matplotlib, например, в виде пустого массива.

Вот мой код на Python:

import osgeo.gdal

dataset = osgeo.gdal.Open("MergedDEM")

gt = dataset.GetGeoTransform()

Answer 1

Вы можете использовать обычный plot_surface метод из матплотлиб. Поскольку ему нужен массив X и Y, он уже нанесен с правильными координатами. Мне всегда трудно создавать красивые 3D-графики, поэтому визуальные аспекты, безусловно, можно улучшить.:)

import gdal
from mpl_toolkits.mplot3d import Axes3D

dem = gdal.Open('gmted_small.tif')
gt  = dem.GetGeoTransform()
dem = dem.ReadAsArray()

fig, ax = plt.subplots(figsize=(16,8), subplot_kw={'projection': '3d'})

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)

X, Y = np.meshgrid(X, Y)

surf = ax.plot_surface(X,Y,dem, rstride=1, cstride=1, cmap=plt.cm.RdYlBu_r, vmin=0, vmax=4000, linewidth=0, antialiased=True)

ax.set_zlim(0, 60000) # to make it stand out less
ax.view_init(60,-105)

fig.colorbar(surf, shrink=0.4, aspect=20)