Python的地形三維可視化Matplotlib和gdal使用實例

更新時間：2017年12月09日 14:42:36 作者：大蝦盧

這篇文章主要介紹了Python的地形三維可視化Matplotlib和gdal使用實例，具有一定借鑒價值，需要的朋友可以了解下。

我是以Python開門的，我還是覺得Python也可以進行地形三維可視化，當然這里需要借助第三方庫，so，我就來介紹：Python一個很重要可視化插件，Matplotlib。

Matplotlib是Python最著名的繪圖庫，它提供了一整套友好的命令，十分適合交互式地進行制圖。而且也可以方便地將它作為繪圖控件，嵌入GUI應用程序中。你會發(fā)現Matplotlib和matlab相似，但是你知道m(xù)atlab強大是很強大，但是安裝包就有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() #最后渲染出你好看的三維圖吧