代码:
# -*- coding: utf-8 -*-
import matplotlib.pylab as plt # 画图模块
import numpy as np
from matplotlib import cm
from matplotlib.colors import LightSource
from osgeo import gdal
raster_path = r'AP_05726_FBS_F0680_RT1.dem.tif'
dataset = gdal.Open(raster_path)
adfGeoTransform = dataset.GetGeoTransform() # 获取投影信息
band = dataset.GetRasterBand(1) # 用gdal去读写你的数据,当然dem只有一个波段
nrows = dataset.RasterXSize
ncols = dataset.RasterYSize # 这两个行就是读取数据的行列数
# 数据的平面四至
lonmin = adfGeoTransform[0]
latmin = adfGeoTransform[3]
lonmax = adfGeoTransform[0] + nrows * adfGeoTransform[1] + ncols * adfGeoTransform[2]
latmax = adfGeoTransform[3] + nrows * adfGeoTransform[4] + ncols * adfGeoTransform[5]
x = np.linspace(lonmin, lonmax, ncols)
y = np.linspace(latmin, latmax, nrows)
# 将数据的x,y,z化作numpy矩阵
lon, lat = np.meshgrid(x, y)
elevation = band.ReadAsArray(0, 0, nrows, ncols)
# 限定一个范围
region = np.s_[10:400, 10:400]
lon, lat, elevation = lon[region], lat[region], elevation[region]
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(12, 10))
ls = LightSource(270, 60) # 设置你可视化数据的色带
rgb = ls.shade(elevation,
cmap=cm.gist_earth, # 设置颜色映射
vert_exag=0.1,
blend_mode='soft')
surf = ax.plot_surface(lon, lat, elevation,
rstride=1, # 指定行的跨度
cstride=1, # 指定列的跨度
facecolors=rgb,
linewidth=0,
antialiased=False, shade=False)
# 设置标题
plt.title("Python-gdal DEM show ", fontsize='large', fontweight='bold', color='#6666FF')
# fig.colorbar(surf, shrink=0.5, aspect=5) # shrink越小,表示colorbar越小
plt.show() # 最后渲染出2.5维图
测试结果: