import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

模拟高程数据（假设数据已经过某种方式插值或生成）

这里我们创建一个简单的40x50网格，并填充随机高程值

x = np.linspace(0, 43.65, 40)
y = np.linspace(0, 58.2, 50)
X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X2 + Y2)) * 100 # 使用一个简单的函数来生成高程数据

绘制三维表面图

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Elevation (m)')
ax.set_title('3D Surface Plot of Elevation Data')

绘制等高线图

plt.subplot(122)
CS = plt.contour(X, Y, Z, colors='k')
plt.clabel(CS, inline=1, fontsize=10)

标注点 A(30,0) 和 B(43,30)

注意：由于X和Y是网格坐标，我们需要找到最接近这些值的索引

idx_a_x = np.argmin(np.abs(x - 30))
idx_a_y = np.argmin(np.abs(y - 0))
idx_b_x = np.argmin(np.abs(x - 43))
idx_b_y = np.argmin(np.abs(y - 30))

plt.plot(x[idx_a_x], y[idx_a_y], 'ro', markersize=5, label='A(30,0)')
plt.plot(x[idx_b_x], y[idx_b_y], 'go', markersize=5, label='B(43,30)')
plt.xlabel('X (km)')
plt.ylabel('Y (km)')
plt.title('Contour Plot of Elevation Data with Points A and B')
plt.legend()

计算地表面积的近似值（忽略地形起伏）

real_area = 43.65 * 58.2
print(f"Actual Surface Area (ignoring elevation changes): {real_area} km^2")

显示图形

plt.tight_layout()
plt.show()

2023310143007