Python的地形三维可视化——简介Matplotlib和gdal

Posted 虾神说D

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Python的地形三维可视化——简介Matplotlib和gdal相关的知识,希望对你有一定的参考价值。

虾神前段时间因为工作关系,去灯塔国享受了一把时差,所以有两周没有更新博客了,幸好有撸文小能手吴道长帮忙,写了这样一篇文章,所以先应急发出来。


题外话:长期不更新还不掉的粉,才是铁杆粉


下面是正文,本文版权属于python代码界后起之秀吴道长所有。


————————我是分割线————————


撸这文,是为了向大虾致敬。大虾曾今用R语言对DEM数据进行三维可视化构建。链接如下:

http://blog.csdn.net/allenlu2008/article/details/49874111

我是以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()  #最后渲染出你好看的三维图吧




最后是虾神的话:GDAL在python里面的安装,属于灰常纠结的事情,运气好(人品好)可能一下就成功,运气不好,可能死活就装不上……如果遇上死活装不上的,可以考虑用MinGW自己进行编译。


如果不想用gdal来实现以上效果,可以考虑用arcpy的RasterToNumpy来把tif文件读成矩阵,也可以实现。


具体还有啥问题的,可以咨询吴道长,吴道长的联系方式,可以给虾神发邮件索取。

以上是关于Python的地形三维可视化——简介Matplotlib和gdal的主要内容,如果未能解决你的问题,请参考以下文章

JAVA三维可视化组件:Matplot 3D for JAVA(V3.0) 一个纯JAVA开发的科学数据可视化组件包 类似 Python 的matplotlib(含示例代码)

基于WebGLCesium技术的三维空间可视化

三维路径规划基于matlab复杂三维地形的无人机路径规划含Matlab源码 2171期

伟景行WebGL为客户创造新价值

Python matplotlib笔记

Cesium入门7 - Adding Terrain - 添加地形