Python的地形三维可视化Matplotlib和gdal使用实例

Z时代
2024-01-10
分类：综合

我是以Python开门的，我还是觉得Python也可以进行地形三维可视化，当然这里需要借助第三方库，so，我就来介绍：Python一个很重要可视化插件，Matplotlib。

Matplotlib是Python最著名的绘图库，它提供了一整套友好的命令，十分适合交互式地进行制图。而且也可以方便地将它作为绘图控件，嵌入GUI应用程序中。你会发现Matplotlib和matlab相似，但是你知道matlab强大是很强大，但是安装包就有7G，一下就让我失去玩弄他的兴趣。

Matplotlib的二维图形非常成熟，可以做各种各样的报表，比如：

好吧，额，....... Matplotlib当然也可以做3维数据。首先准备了一份样例，珠穆朗姆峰周围地势的dem数据，我在arcsense进行拉伸显示,初步看看效果：

讲这份数据，输入编写好的相关路径，然后经过简单的三维渲染：

Matplotlib只是一个可视化包，不能直接识别栅格数据，所以这边还用到另外一个包---gdal，开源栅格空间数据转换库，他的本体竟然还是c++，也开放了一个Python框架，这算是向Python投诚嘛？

好吧，开源一下用Matplotlib和gdal包进行地形三维可视化的小脚本：


# -*- coding: gbk -*-
from mpl_toolkits.mplot3d import Axes3D
from Matplotlib import cbook
from Matplotlib import cm
from Matplotlib.colors import LightSource
import Matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal
gdal.AllRegister()
filePath = u"E:\桌面\dem实验/Himalaya.tif" #输入你的dem数据
dataset = gdal.Open(filePath)
adfGeoTransform = dataset.GetGeoTransform()
band = dataset.GetRasterBand(1)  #用gdal去读写你的数据，当然dem只有一个波段
nrows = dataset.RasterXSize 
ncols  = dataset.RasterYSize  #这两个行就是读取数据的行列数
Xmin = adfGeoTransform[0] #你的数据的平面四至
Ymin = adfGeoTransform[3]
Xmax = adfGeoTransform[0] + nrows * adfGeoTransform[1] + ncols * adfGeoTransform[2]
Ymax = adfGeoTransform[3] + nrows * adfGeoTransform[4] + ncols * adfGeoTransform[5]
x = np.linspace(Xmin,Xmax, ncols)
y = np.linspace(Ymin,Ymax, nrows)
X,Y = np.meshgrid(x, y)
Z = band.ReadAsArray(0, 0,nrows, ncols) #这一段就是讲数据的x，y，z化作numpy矩阵
region = np.s_[10:400,10:400]
X, Y, Z = X[region], Y[region],Z[region]
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(12,10))
ls = LightSource(270, 20)  #设置你可视化数据的色带
rgb = ls.shade(Z, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=rgb,
            linewidth=0, antialiased=False, shade=False)
plt.show() #最后渲染出你好看的三维图吧