首先请大家读一下面这篇文章了解什么是Gdal

http://blog.csdn.net/grllery/article/details/77822595

剩下的我要公布绘制富士山的代码了，虽然基本copy虾神的路子，我加入一些注释方便理解

# -*- coding: utf-8 -*-

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 = "ceshi.tif" #输入你的dem数据(我放在工作目录，你要记得修改你的哦O^-^O)

dataset = gdal.Open(filePath)

adfGeoTransform = dataset.GetGeoTransform()

band = dataset.GetRasterBand(1) #用gdal去读写你的数据，当然dem只有一个波段

ncols = dataset.RasterXSize #图像的宽度(X方向上的像素个数) 数据的列数 (这里就是Gdal中的格子和矩阵的不同了，这个疑问只要你好好读上面的那篇文章，不难理解哈哈O^-^O)

nrows = dataset.RasterYSize#图像的宽度(Y方向上的像素个数) 数据的行数

Xmin = adfGeoTransform[0] #你的数据的平面四至

Ymin = adfGeoTransform[3]

Xmax = adfGeoTransform[0] + nrows * adfGeoTransform[1] + ncols * adfGeoTransform[2]

Ymax = adfGeoTransform[3] + nrows * adfGeoTransform[4] + ncols * adfGeoTransform[5]#(这几个参数也是在那篇文章中介绍了O^-^O)

x = np.linspace(Xmin,Xmax, ncols)#地理x坐标 数组y坐标

y = np.linspace(Ymin,Ymax, nrows)#地理y坐标 数组x坐标

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

Z = band.ReadAsArray(0, 0,ncols, nrows) #这一段就是讲数据的x，y，z化作numpy矩阵(这里Z读取完后是一个Y方向的像素个数*X方向上的像素个数的一个矩阵O^-^O)

region = np.s_[10:400,10:400] #这家伙就等同于一个切片命令(是的没错就这货slice(start,stop,step)         23333333333333)

X, Y, Z = X[region], Y[region],Z[region]#数组转置和轴对换：数组不仅有transpose方法，还有一个特殊的T属性比如Z[region].T

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

ps：我用的python是Anaconda集成的，避免了安装Matplot         (233333333333)

Anaconda中安装Gdal，安装方法灰常简单 WIN+R cmd进入命令窗口然后输入python回车，pip install gdal，让他自己安去吧，静待successful

我的联系方式 Email：mafengkai@yahoo.com  不对的地方还请大家提醒了  BYEBYE！