[python]gdal安装后测试代码

# -*- 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维图