2.6 某区域43.65km*58.2km的高程数据，画出该区域的三维表面图和等高线图，在A（30，0）和B（43，30）（单位：km）点处建立了两个基地，在等高线图上标注出这两个点，并求该区域地表面积的近似值

点击查看代码

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(X**2 + Y**2)) * 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()
 
print("学号：3004")