Отображение поверхности DEM с привязкой в ​​3D matplotlib

Question

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

Вот мой питон код до сих пор:

import osgeo.gdal 

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

gt = dataset.GetGeoTransform() 

Answer 1

Вы можете использовать обычный метод plot_surface из Matplotlib. Потому что ему нужен массив 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)