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

"""
@File    : 02_读取dem影像2.5维展示.py
@Author  : fungis@163.com
@Time    : 2022/09/17 23:32
@notice  :
"""

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'./data-use/tif/AP_05726_FBS_F0680_RT1.dem.tif'
dataset = gdal.Open(raster_path)
data_transform = dataset.GetGeoTransform()  # 获取投影信息
band = dataset.GetRasterBand(1)  # 用gdal去读写你的数据，当然dem只有一个波段

nrows = dataset.RasterXSize
ncols = dataset.RasterYSize  # 这两个行就是读取数据的行列数

# 数据的平面四至
lon_min = data_transform[0]
lat_min = data_transform[3]
lon_max = data_transform[0] + nrows * data_transform[1] + ncols * data_transform[2]
lat_max = data_transform[3] + nrows * data_transform[4] + ncols * data_transform[5]

x = np.linspace(lon_min, lon_max, ncols)
y = np.linspace(lat_min, lat_max, 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,  # cmap是颜色映射表
               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)
''' 等高线图绘制
ax.contour(lon, lat, elevation, zdir='elevation', offset=-2, cmap=plt.get_cmap('rainbow')) 
# 设置z轴的维度，x,y类似
ax.set_zlim(-2, 2)
'''

# 设置标题
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维图