python实现dem输出三维模型_python实现DEM数据的阴影生成的方法

要用 Python 实现 DEM 输出三维模型可以使用一些开源的库和工具，比如：

1. GDAL：可以用于读取 DEM 数据，生成高程网格等。
2. PyVista：可以用于可视化 DEM 数据，并生成三维模型。
3. Blender：可以通过 Python 脚本生成三维模型。

以下是使用 PyVista 库生成三维模型的示例代码：

import pyvista as pv
import rasterio as rio

# 读取 DEM 文件
with rio.open('dem.tif') as src:
    z = src.read(1)
    transform = src.transform

# 创建 PyVista 网格
grid = pv.UnstructuredGrid(pv.vtk.vtkXMLUnstructuredGridReader().GetOutput())
grid.points = pv.vtk.vtkPoints()
grid.cells = pv.vtk.vtkCellArray()

# 将 DEM 中每个像素转换为网格节点
for i in range(z.shape[0]):
    for j in range(z.shape[1]):
        grid.points.InsertNextPoint(transform * (j, i, z[i, j]))

# 将节点连接成三角形网格
for i in range(z.shape[0] - 1):
    for j in range(z.shape[1] - 1):
        cell = pv.vtk.vtkTriangle()
        cell.GetPointIds().SetId(0, i * z.shape[1] + j)
        cell.GetPointIds().SetId(1, i * z.shape[1] + j + 1)
        cell.GetPointIds().SetId(2, (i + 1) * z.shape[1] + j)
        grid.cells.InsertNextCell(cell)

# 可视化网格
plotter = pv.Plotter()
plotter.add_mesh(grid, cmap='terrain')
plotter.show()

要用 Python 实现 DEM 数据的阴影生成，可以使用以下步骤：

1. 计算光照方向。可以根据 DEM 数据的法向量和太阳高度角、方位角等参数计算。
2. 计算每个像素的阴影值。可以使用遮蔽法或者投影法等方法。
3. 可视化 DEM 数据，并根据阴影值着色。可以使用 PyVista 或者 Matplotlib 等库进行可视化。

以下是使用遮蔽法计算阴影值的示例代码：

import numpy as np
import rasterio as rio

# 读取 DEM 文件
with rio.open('dem.tif') as src:
    z = src.read(1)
    transform = src.transform

# 计算 DEM 数据的法向量
dx, dy = np.gradient(z, transform[0], transform[4])
dz = np.ones_like(z)
nx, ny, nz = np.cross(dx, dy, dz, axis=0, normalize=True)

# 计算光照方向
theta = np.radians(30)  # 太阳高度角
phi = np.radians(-45)  # 太阳方位角
lx = np.cos(theta) * np.cos(phi)
ly = np.cos(theta) * np.sin(phi)
lz = np.sin(theta)

# 计算每个像素的阴影值
shadow = np.zeros_like(z, dtype=bool)
for i in range(z.shape[0]):
    for j in range(z.shape[1]):
        p = (j, i, z[i, j])
        q = (p[0] - lx, p[1] - ly, p[2] - lz)
        r = (p[0] + lx, p[1] + ly, p[2] + lz)
        if np.dot(nz[i, j], q - p) > 0 and np.dot(nz[i, j], r - p) > 0:
            shadow[i, j] = True

# 可视化 DEM 数据，并根据阴影值着色
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

fig, ax = plt.subplots()
ls = LightSource(azdeg=315, altdeg=45)
rgb = ls.shade(z, cmap='terrain', vert_exag=10, blend_mode='soft', dx=transform[0], dy=-transform[4], fraction=1, shade_alpha=0.8)
ax.imshow(rgb)
ax.imshow(shadow, cmap='gray', alpha=0.5)
plt.show()