首先请大家读一下面这篇文章了解什么是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() #最后渲染出你好看的三维图吧（）
ps：我用的python是Anaconda集成的，避免了安装Matplot （）
Anaconda中安装Gdal，安装方法灰常简单 WIN+R cmd进入命令窗口然后输入python回车，pip install gdal，让他自己安去吧，静待successful
我的联系方式 Email：mafengkai@yahoo.com 不对的地方还请大家提醒了 BYEBYE！