import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib.font_manager as fm
from scipy.interpolate import griddata, Rbf
from scipy.ndimage import gaussian_filter
import trimesh
import os
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
class SlopeTemperatureAnalysis:
def __init__(self, temp_data_file, coord_data_file, stl_file,
xy_slice_z=None, xz_slice_y=None, yz_slice_x=None):
"""
初始化边坡温度分析系统
参数:
temp_data_file: 测点-时间-温度变化表Excel文件路径
coord_data_file: 测点坐标表Excel文件路径
stl_file: 边坡三维STL文件路径
xy_slice_z: XY平面切片的Z坐标(mm)，如果为None则使用平均高度
xz_slice_y: XZ平面切片的Y坐标(mm)，如果为None则使用中心位置
yz_slice_x: YZ平面切片的X坐标(mm)，如果为None则使用中心位置
"""
# 设置中文字体
self.set_chinese_font()
# 设置工作目录
self.base_dir = r'D:\桌面\温度数据\三维温度图\工况1'
# 创建输出文件夹
self.output_dir = os.path.join(self.base_dir, '分析结果')
os.makedirs(self.output_dir, exist_ok=True)
# 切片位置设置
self.xy_slice_z = xy_slice_z
self.xz_slice_y = xz_slice_y
self.yz_slice_x = yz_slice_x
# 构建完整文件路径
temp_data_path = os.path.join(self.base_dir, temp_data_file)
coord_data_path = os.path.join(self.base_dir, coord_data_file)
stl_path = os.path.join(self.base_dir, stl_file)
# 检查文件是否存在
self.check_files_exist(temp_data_path, coord_data_path, stl_path)
# 加载数据
self.load_data(temp_data_path, coord_data_path)
# 加载边坡模型
self.load_slope_model(stl_path)
# 预处理数据
self.preprocess_data()
print("边坡温度分析系统初始化完成")
print(f"结果将保存到: {self.output_dir}")
print(f"切片位置设置 - XY平面Z坐标: {self.xy_slice_z}mm, XZ平面Y坐标: {self.xz_slice_y}mm, YZ平面X坐标: {self.yz_slice_x}mm")
def check_files_exist(self, temp_path, coord_path, stl_path):
"""检查必要的文件是否存在"""
if not os.path.exists(temp_path):
print(f"警告: 温度数据文件不存在: {temp_path}")
if not os.path.exists(coord_path):
print(f"警告: 坐标数据文件不存在: {coord_path}")
if not os.path.exists(stl_path):
print(f"警告: STL模型文件不存在: {stl_path}")
def set_chinese_font(self):
"""设置中文字体支持"""
try:
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
except:
plt.rcParams['font.sans-serif'] = ['DejaVu Sans']
def load_data(self, temp_data_file, coord_data_file):
"""
从Excel文件加载温度数据和坐标数据
"""
print("正在从Excel文件加载数据...")
try:
# 加载温度数据
self.temp_df = pd.read_excel(temp_data_file)
print(f"温度数据记录数: {len(self.temp_df)}")
print(f"温度数据列名: {list(self.temp_df.columns)}")
# 加载坐标数据
self.coord_df = pd.read_excel(coord_data_file)
print(f"测点数量: {len(self.coord_df)}")
print(f"坐标数据列名: {list(self.coord_df.columns)}")
except Exception as e:
print(f"Excel文件读取失败: {e}")
print("尝试使用CSV格式...")
# 如果Excel读取失败，尝试CSV格式
self.temp_df = pd.read_csv(temp_data_file)
self.coord_df = pd.read_csv(coord_data_file)
# 检查并标准化列名
self.standardize_column_names()
# 处理重复列名
self.handle_duplicate_columns()
# 合并数据
self.merged_df = pd.merge(self.temp_df, self.coord_df, on='测点编号')
print(f"合并后数据记录数: {len(self.merged_df)}")
# 转换时间格式
if '时间' in self.merged_df.columns:
self.merged_df['时间'] = pd.to_datetime(self.merged_df['时间'])
elif '日期' in self.merged_df.columns:
self.merged_df['时间'] = pd.to_datetime(self.merged_df['日期'])
print("使用'日期'列作为时间列")
def standardize_column_names(self):
"""标准化列名"""
# 温度数据列名映射
temp_column_mapping = {}
for col in self.temp_df.columns:
col_lower = str(col).lower()
if any(keyword in col_lower for keyword in ['测点', '点号', 'point']):
temp_column_mapping[col] = '测点编号'
elif any(keyword in col_lower for keyword in ['时间', '日期', 'time', 'date']):
temp_column_mapping[col] = '时间'
elif any(keyword in col_lower for keyword in ['温度', 'temp']):
temp_column_mapping[col] = '温度'
if temp_column_mapping:
self.temp_df = self.temp_df.rename(columns=temp_column_mapping)
print(f"温度数据列名标准化: {temp_column_mapping}")
# 坐标数据列名映射
coord_column_mapping = {}
for col in self.coord_df.columns:
col_lower = str(col).lower()
if any(keyword in col_lower for keyword in ['测点', '点号', 'point']):
coord_column_mapping[col] = '测点编号'
elif any(keyword in col_lower for keyword in ['x', '坐标x', 'coordx']):
coord_column_mapping[col] = 'X坐标'
elif any(keyword in col_lower for keyword in ['y', '坐标y', 'coordy']):
coord_column_mapping[col] = 'Y坐标'
elif any(keyword in col_lower for keyword in ['高程', '海拔', 'z', '坐标z', 'coordz', 'elevation']):
coord_column_mapping[col] = '高程'
if coord_column_mapping:
self.coord_df = self.coord_df.rename(columns=coord_column_mapping)
print(f"坐标数据列名标准化: {coord_column_mapping}")
# 检查必要的列是否存在
required_temp_cols = ['测点编号', '时间', '温度']
required_coord_cols = ['测点编号', 'X坐标', 'Y坐标', '高程']
missing_temp = [col for col in required_temp_cols if col not in self.temp_df.columns]
missing_coord = [col for col in required_coord_cols if col not in self.coord_df.columns]
if missing_temp:
print(f"警告: 温度数据缺少必要的列: {missing_temp}")
if missing_coord:
print(f"警告: 坐标数据缺少必要的列: {missing_coord}")
def handle_duplicate_columns(self):
"""处理重复的列名"""
# 检查温度数据中的重复列名
temp_duplicates = self.temp_df.columns.duplicated()
if temp_duplicates.any():
print(f"温度数据中有重复列名: {list(self.temp_df.columns[temp_duplicates])}")
# 删除重复列，保留第一个
self.temp_df = self.temp_df.loc[:, ~temp_duplicates]
print("已删除温度数据中的重复列")
# 检查坐标数据中的重复列名
coord_duplicates = self.coord_df.columns.duplicated()
if coord_duplicates.any():
print(f"坐标数据中有重复列名: {list(self.coord_df.columns[coord_duplicates])}")
# 删除重复列，保留第一个
self.coord_df = self.coord_df.loc[:, ~coord_duplicates]
print("已删除坐标数据中的重复列")
# 打印最终列名
print(f"温度数据最终列名: {list(self.temp_df.columns)}")
print(f"坐标数据最终列名: {list(self.coord_df.columns)}")
def load_slope_model(self, stl_file):
"""
从STL文件加载边坡三维模型
"""
print("正在从STL文件加载边坡三维模型...")
try:
# 使用trimesh加载STL文件
self.slope_mesh = trimesh.load_mesh(stl_file)
print(f"STL文件加载成功:")
print(f"顶点数: {len(self.slope_mesh.vertices)}")
print(f"面片数: {len(self.slope_mesh.faces)}")
print(f"模型边界: {self.slope_mesh.bounds}")
# 获取模型的几何中心
self.mesh_center = self.slope_mesh.centroid
print(f"模型中心: {self.mesh_center}")
# 检查模型是否是水密的（封闭的）
self.is_watertight = self.slope_mesh.is_watertight
print(f"模型是否水密: {self.is_watertight}")
except Exception as e:
print(f"STL文件读取失败: {e}")
print("创建基于测点坐标的简化边坡模型...")
self.create_simplified_slope_model()
def create_simplified_slope_model(self):
"""创建基于测点坐标的简化边坡模型"""
print("创建基于测点坐标的简化边坡模型...")
# 使用测点坐标创建凸包作为简化模型
points = self.coord_df[['X坐标', 'Y坐标', '高程']].values
self.slope_mesh = trimesh.convex.convex_hull(points)
self.is_watertight = self.slope_mesh.is_watertight
print(f"简化模型边界: {self.slope_mesh.bounds}")
def is_point_in_mesh(self, points, tolerance=0.1):
"""
判断点是否在网格模型内部或边界附近
参数:
points: (n, 3) 形状的点坐标数组
tolerance: 边界容差(mm)，用于保留边界附近的点
返回:
inside: (n,) 形状的布尔数组，True表示点在模型内部或边界附近
"""
if not hasattr(self, 'slope_mesh') or self.slope_mesh is None:
# 如果没有模型，认为所有点都在范围内
return np.ones(len(points), dtype=bool)
try:
# 使用射线法判断点是否在模型内部
if self.is_watertight:
inside = self.slope_mesh.contains(points)
# 对于不在内部的点，检查是否在边界附近
if not np.all(inside):
# 计算到模型表面的距离
distances = self.slope_mesh.nearest.on_surface(points[~inside])[1]
# 如果距离小于容差，则认为在边界上
on_boundary = distances < tolerance
inside[~inside] = on_boundary
else:
# 如果模型不是水密的，使用近似方法：判断点是否在模型边界框内
bounds = self.slope_mesh.bounds
inside = np.all((points >= bounds[0]) & (points <= bounds[1]), axis=1)
return inside
except Exception as e:
print(f"点包含性检测失败: {e}")
# 失败时返回所有点都在范围内
return np.ones(len(points), dtype=bool)
def preprocess_data(self):
"""
数据预处理
"""
print("正在预处理数据...")
# 提取唯一时间点（包括日期和时间）
if '时间' in self.merged_df.columns:
# 提取唯一的时间点
self.time_points = sorted(self.merged_df['时间'].unique())
print(f"数据覆盖时间点: {len(self.time_points)}")
# 同时提取日期用于分组
self.merged_df['日期'] = self.merged_df['时间'].dt.date
self.dates = sorted(self.merged_df['日期'].unique())
print(f"数据覆盖天数: {len(self.dates)}")
else:
print("警告: 未找到时间列，使用默认日期")
self.time_points = [datetime.now()]
self.dates = [datetime.now().date()]
# 提取测点信息
self.monitoring_points = self.merged_df[['测点编号', 'X坐标', 'Y坐标', '高程']].drop_duplicates()
print(f"监测点数量: {len(self.monitoring_points)}")
# 创建插值网格（只在模型范围内）
self.create_interpolation_grid()
def create_interpolation_grid(self):
"""
创建温度插值网格 - 只在STL模型范围内创建网格点
"""
# 如果有STL模型，使用模型边界
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
min_x, max_x = bounds[0][0], bounds[1][0]
min_y, max_y = bounds[0][1], bounds[1][1]
min_z, max_z = bounds[0][2], bounds[1][2]
else:
# 否则使用测点坐标范围
min_x, max_x = self.coord_df['X坐标'].min(), self.coord_df['X坐标'].max()
min_y, max_y = self.coord_df['Y坐标'].min(), self.coord_df['Y坐标'].max()
min_z, max_z = self.coord_df['高程'].min(), self.coord_df['高程'].max()
# 创建密集的初始网格
grid_points = 50
# 将网格数组保存为实例属性
self.x_grid = np.linspace(min_x, max_x, grid_points)
self.y_grid = np.linspace(min_y, max_y, grid_points)
self.z_grid = np.linspace(min_z, max_z, grid_points // 2)
# 创建完整的网格
self.X_grid, self.Y_grid, self.Z_grid = np.meshgrid(
self.x_grid, self.y_grid, self.z_grid, indexing='ij'
)
# 将网格点展平以进行包含性检测
grid_points_flat = np.vstack([
self.X_grid.ravel(),
self.Y_grid.ravel(),
self.Z_grid.ravel()
]).T
# 检测哪些点在模型内部或边界附近
print("正在检测网格点在模型内部的情况...")
self.inside_mask = self.is_point_in_mesh(grid_points_flat, tolerance=0.5)
inside_ratio = np.sum(self.inside_mask) / len(self.inside_mask)
print(f"模型内部网格点比例: {inside_ratio:.2%}")
print(f"有效网格点数: {np.sum(self.inside_mask)}")
print(f"插值网格尺寸: {self.X_grid.shape}")
print(f"网格范围: X({min_x:.1f}-{max_x:.1f}mm), Y({min_y:.1f}-{max_y:.1f}mm), Z({min_z:.1f}-{max_z:.1f}mm)")
def format_time_point(self, time_point):
"""
格式化时间点，处理numpy.datetime64和datetime对象
格式: 月-日 小时 (省略年份，精确到小时)
"""
if isinstance(time_point, np.datetime64):
# 将numpy.datetime64转换为pandas.Timestamp
time_point = pd.Timestamp(time_point)
if hasattr(time_point, 'strftime'):
# 格式化为 "月-日 小时" 格式，例如 "07-29 15"
return time_point.strftime("%m-%d %H")
else:
return str(time_point)
def interpolate_temperature_3d(self, time_point):
"""
三维温度场插值 - 只在模型范围内插值
参数:
time_point: 时间点
返回:
三维温度场（模型外部的点为NaN）
"""
# 获取指定时间点的数据
time_data = self.merged_df[self.merged_df['时间'] == time_point]
if len(time_data) == 0:
print(f"警告: 时间点 {time_point} 无数据")
return None
# 提取测点坐标和温度
points = time_data[['X坐标', 'Y坐标', '高程']].values
temperatures = time_data['温度'].values
print(f"插值使用测点数: {len(points)}")
# 使用径向基函数插值获得更平滑的结果
try:
# 尝试使用RBF插值
rbf = Rbf(points[:, 0], points[:, 1], points[:, 2], temperatures,
function='linear', smooth=0.1)
# 在网格上计算温度
grid_temps = rbf(self.X_grid, self.Y_grid, self.Z_grid)
# 应用高斯滤波平滑结果
grid_temps = gaussian_filter(grid_temps, sigma=0.5)
except Exception as e:
print(f"RBF插值失败: {e}，使用线性插值")
# 如果RBF失败，使用线性插值
grid_temps = griddata(
points, temperatures,
(self.X_grid, self.Y_grid, self.Z_grid),
method='linear', fill_value=np.nan
)
# 将模型外部的温度设为NaN
grid_temps_flat = grid_temps.ravel()
grid_temps_flat[~self.inside_mask] = np.nan
grid_temps = grid_temps_flat.reshape(grid_temps.shape)
return grid_temps
def create_time_temperature_contour(self, time_point):
"""
创建每个时间点的温度场等值线云图
"""
time_str = self.format_time_point(time_point)
print(f"正在生成 {time_str} 的温度场等值线云图...")
# 获取温度场数据
temp_field = self.interpolate_temperature_3d(time_point)
if temp_field is None:
return
# 创建可视化
fig = plt.figure(figsize=(18, 12))
# 1. 三维等值面图
ax1 = fig.add_subplot(231, projection='3d')
self.plot_3d_isosurface(ax1, temp_field, time_point)
# 2. 三维切片图
ax2 = fig.add_subplot(232, projection='3d')
self.plot_3d_slices(ax2, temp_field, time_point)
# 3. XY平面温度等值线
ax3 = fig.add_subplot(233)
self.plot_xy_contour(ax3, temp_field, time_point)
# 4. XZ平面温度等值线
ax4 = fig.add_subplot(234)
self.plot_xz_contour(ax4, temp_field, time_point)
# 5. YZ平面温度等值线
ax5 = fig.add_subplot(235)
self.plot_yz_contour(ax5, temp_field, time_point)
# 6. 温度分布直方图
ax6 = fig.add_subplot(236)
self.plot_temperature_histogram(ax6, temp_field, time_point)
plt.tight_layout()
# 保存图片
filename = f"温度场等值线_{self.format_time_point(time_point).replace(':', '').replace(' ', '_')}.png"
filepath = os.path.join(self.output_dir, filename)
plt.savefig(filepath, dpi=300, bbox_inches='tight')
plt.close()
print(f"已保存: {filename}")
def plot_3d_isosurface(self, ax, temp_field, time_point):
"""
绘制三维等值面图 - 只在模型范围内绘制，网格边界完全透明
"""
try:
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
# 创建掩码，只处理模型内部的点
valid_mask = ~np.isnan(temp_field)
if np.sum(valid_mask) == 0:
print("警告: 没有有效的温度数据用于等值面生成")
return
# 计算等值面水平
valid_temps = temp_field[valid_mask]
temp_min = np.nanmin(valid_temps)
temp_max = np.nanmax(valid_temps)
levels = np.linspace(temp_min, temp_max, 6)
# 创建用于marching cubes的数组（将NaN设为很小的值）
temp_field_mc = np.copy(temp_field)
temp_field_mc[np.isnan(temp_field_mc)] = temp_min - 10
# 绘制等值面
for i, level in enumerate(levels[1:-1]): # 跳过最低和最高值
try:
# 计算网格间距
dx = self.x_grid[1] - self.x_grid[0] if len(self.x_grid) > 1 else 1.0
dy = self.y_grid[1] - self.y_grid[0] if len(self.y_grid) > 1 else 1.0
dz = self.z_grid[1] - self.z_grid[0] if len(self.z_grid) > 1 else 1.0
# 使用 marching cubes 算法提取等值面
verts, faces, _, _ = measure.marching_cubes(
temp_field_mc, level=level, spacing=(dx, dy, dz)
)
# 转换坐标到实际位置
verts[:, 0] = verts[:, 0] + self.x_grid[0]
verts[:, 1] = verts[:, 1] + self.y_grid[0]
verts[:, 2] = verts[:, 2] + self.z_grid[0]
# 创建多边形集合 - 设置网格边界完全透明
mesh = Poly3DCollection(verts[faces], alpha=0.3, linewidth=0)
# 根据温度值设置颜色
normalized_temp = (level - temp_min) / (temp_max - temp_min)
color = plt.cm.jet(normalized_temp)
mesh.set_facecolor(color)
mesh.set_edgecolor('none') # 设置边缘颜色为完全透明
ax.add_collection3d(mesh)
except Exception as e:
print(f"等值面 {level:.1f}°C 生成失败: {e}")
continue
except ImportError:
print("skimage不可用，跳过等值面生成")
except Exception as e:
print(f"等值面生成失败: {e}")
# 绘制监测点
time_data = self.merged_df[self.merged_df['时间'] == time_point]
scatter = ax.scatter(time_data['X坐标'], time_data['Y坐标'], time_data['高程'],
c=time_data['温度'], cmap='jet', s=50,
edgecolors='black', linewidth=1, alpha=0.8)
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
ax.set_xlim(bounds[0][0], bounds[1][0])
ax.set_ylim(bounds[0][1], bounds[1][1])
ax.set_zlim(bounds[0][2], bounds[1][2])
ax.set_xlabel('X坐标 (mm)', fontsize=10)
ax.set_ylabel('Y坐标 (mm)', fontsize=10)
ax.set_zlabel('高程 (mm)', fontsize=10)
ax.set_title(f'三维温度等值面\n{self.format_time_point(time_point)}', fontsize=12)
# 添加颜色条
plt.colorbar(scatter, ax=ax, shrink=0.6, label='温度 (°C)')
def plot_3d_slices(self, ax, temp_field, time_point):
"""
绘制三维切片图 - 只在模型范围内绘制，网格边界完全透明
"""
# 选择几个切片位置
x_slice_idx = len(self.x_grid) // 2
y_slice_idx = len(self.y_grid) // 2
z_slice_idx = len(self.z_grid) // 3
# 创建颜色归一化
valid_temps = temp_field[~np.isnan(temp_field)]
if len(valid_temps) == 0:
print("警告: 没有有效的温度数据")
return
temp_min = np.min(valid_temps)
temp_max = np.max(valid_temps)
norm = plt.Normalize(temp_min, temp_max)
# X方向切片
xx, yy = np.meshgrid(self.y_grid, self.z_grid)
zz = np.ones_like(xx) * self.x_grid[x_slice_idx]
slice_temp_x = temp_field[x_slice_idx, :, :].T
# 只绘制有效数据
valid_x = ~np.isnan(slice_temp_x)
if np.any(valid_x):
colors_x = plt.cm.jet(norm(slice_temp_x))
# 将无效数据设为透明
colors_x[~valid_x] = [0, 0, 0, 0]
surf_x = ax.plot_surface(zz, xx, yy, facecolors=colors_x,
alpha=0.6, rstride=2, cstride=2, linewidth=0, edgecolor='none')
# Y方向切片
xx, zz = np.meshgrid(self.x_grid, self.z_grid)
yy = np.ones_like(xx) * self.y_grid[y_slice_idx]
slice_temp_y = temp_field[:, y_slice_idx, :].T
valid_y = ~np.isnan(slice_temp_y)
if np.any(valid_y):
colors_y = plt.cm.jet(norm(slice_temp_y))
colors_y[~valid_y] = [0, 0, 0, 0]
surf_y = ax.plot_surface(xx, yy, zz, facecolors=colors_y,
alpha=0.6, rstride=2, cstride=2, linewidth=0, edgecolor='none')
# Z方向切片（地表）
xx, yy = np.meshgrid(self.x_grid, self.y_grid)
zz = np.ones_like(xx) * self.z_grid[z_slice_idx]
slice_temp_z = temp_field[:, :, z_slice_idx].T
valid_z = ~np.isnan(slice_temp_z)
if np.any(valid_z):
colors_z = plt.cm.jet(norm(slice_temp_z))
colors_z[~valid_z] = [0, 0, 0, 0]
surf_z = ax.plot_surface(xx, yy, zz, facecolors=colors_z,
alpha=0.7, rstride=2, cstride=2, linewidth=0, edgecolor='none')
# 绘制监测点
time_data = self.merged_df[self.merged_df['时间'] == time_point]
scatter = ax.scatter(time_data['X坐标'], time_data['Y坐标'], time_data['高程'],
c=time_data['温度'], cmap='jet', s=50,
edgecolors='white', linewidth=1)
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
ax.set_xlim(bounds[0][0], bounds[1][0])
ax.set_ylim(bounds[0][1], bounds[1][1])
ax.set_zlim(bounds[0][2], bounds[1][2])
ax.set_xlabel('X坐标 (mm)', fontsize=10)
ax.set_ylabel('Y坐标 (mm)', fontsize=10)
ax.set_zlabel('高程 (mm)', fontsize=10)
ax.set_title(f'三维温度切片\n{self.format_time_point(time_point)}', fontsize=12)
# 添加颜色条
plt.colorbar(scatter, ax=ax, shrink=0.6, label='温度 (°C)')
def plot_xy_contour(self, ax, temp_field, time_point):
"""
绘制XY平面温度等值线 - 只在模型范围内绘制
"""
# 确定切片位置
if self.xy_slice_z is not None:
# 使用用户指定的Z坐标
z_idx = np.argmin(np.abs(self.z_grid - self.xy_slice_z))
slice_z = self.xy_slice_z
else:
# 使用平均高度
z_idx = len(self.z_grid) // 2
slice_z = self.z_grid[z_idx]
slice_temp = temp_field[:, :, z_idx]
# 创建等值线图（只绘制有效数据）
valid_mask = ~np.isnan(slice_temp)
if np.any(valid_mask):
contour = ax.contourf(self.x_grid, self.y_grid, slice_temp.T,
levels=20, cmap='jet', alpha=0.8)
# 添加等值线
contour_lines = ax.contour(self.x_grid, self.y_grid, slice_temp.T,
levels=10, colors='black', linewidths=0.5, alpha=0.7)
ax.clabel(contour_lines, inline=True, fontsize=8, fmt='%.1f°C')
# 添加监测点
time_data = self.merged_df[self.merged_df['时间'] == time_point]
scatter = ax.scatter(time_data['X坐标'], time_data['Y坐标'],
c=time_data['温度'], cmap='jet', s=50,
edgecolors='white', linewidth=1)
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
ax.set_xlim(bounds[0][0], bounds[1][0])
ax.set_ylim(bounds[0][1], bounds[1][1])
ax.set_xlabel('X坐标 (mm)', fontsize=10)
ax.set_ylabel('Y坐标 (mm)', fontsize=10)
ax.set_title(f'XY平面温度等值线\n高程: {slice_z:.1f}mm\n{self.format_time_point(time_point)}', fontsize=12)
ax.grid(True, alpha=0.3)
if np.any(valid_mask):
plt.colorbar(contour, ax=ax, label='温度 (°C)')
def plot_xz_contour(self, ax, temp_field, time_point):
"""
绘制XZ平面温度等值线 - 只在模型范围内绘制
"""
# 确定切片位置
if self.xz_slice_y is not None:
# 使用用户指定的Y坐标
y_idx = np.argmin(np.abs(self.y_grid - self.xz_slice_y))
slice_y = self.xz_slice_y
else:
# 使用中心位置
y_idx = len(self.y_grid) // 2
slice_y = self.y_grid[y_idx]
slice_temp = temp_field[:, y_idx, :]
# 创建等值线图（只绘制有效数据）
valid_mask = ~np.isnan(slice_temp)
if np.any(valid_mask):
contour = ax.contourf(self.x_grid, self.z_grid, slice_temp.T,
levels=20, cmap='jet', alpha=0.8)
# 添加等值线
contour_lines = ax.contour(self.x_grid, self.z_grid, slice_temp.T,
levels=10, colors='black', linewidths=0.5, alpha=0.7)
ax.clabel(contour_lines, inline=True, fontsize=8, fmt='%.1f°C')
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
ax.set_xlim(bounds[0][0], bounds[1][0])
ax.set_ylim(bounds[0][2], bounds[1][2])
ax.set_xlabel('X坐标 (mm)', fontsize=10)
ax.set_ylabel('高程 (mm)', fontsize=10)
ax.set_title(f'XZ平面温度等值线\nY: {slice_y:.1f}mm\n{self.format_time_point(time_point)}', fontsize=12)
ax.grid(True, alpha=0.3)
if np.any(valid_mask):
plt.colorbar(contour, ax=ax, label='温度 (°C)')
def plot_yz_contour(self, ax, temp_field, time_point):
"""
绘制YZ平面温度等值线 - 只在模型范围内绘制
"""
# 确定切片位置
if self.yz_slice_x is not None:
# 使用用户指定的X坐标
x_idx = np.argmin(np.abs(self.x_grid - self.yz_slice_x))
slice_x = self.yz_slice_x
else:
# 使用中心位置
x_idx = len(self.x_grid) // 2
slice_x = self.x_grid[x_idx]
slice_temp = temp_field[x_idx, :, :]
# 创建等值线图（只绘制有效数据）
valid_mask = ~np.isnan(slice_temp)
if np.any(valid_mask):
contour = ax.contourf(self.y_grid, self.z_grid, slice_temp.T,
levels=20, cmap='jet', alpha=0.8)
# 添加等值线
contour_lines = ax.contour(self.y_grid, self.z_grid, slice_temp.T,
levels=10, colors='black', linewidths=0.5, alpha=0.7)
ax.clabel(contour_lines, inline=True, fontsize=8, fmt='%.1f°C')
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
ax.set_xlim(bounds[0][1], bounds[1][1])
ax.set_ylim(bounds[0][2], bounds[1][2])
ax.set_xlabel('Y坐标 (mm)', fontsize=10)
ax.set_ylabel('高程 (mm)', fontsize=10)
ax.set_title(f'YZ平面温度等值线\nX: {slice_x:.1f}mm\n{self.format_time_point(time_point)}', fontsize=12)
ax.grid(True, alpha=0.3)
if np.any(valid_mask):
plt.colorbar(contour, ax=ax, label='温度 (°C)')
def plot_temperature_histogram(self, ax, temp_field, time_point):
"""
绘制温度分布直方图
"""
# 获取温度数据
time_data = self.merged_df[self.merged_df['时间'] == time_point]
measured_temps = time_data['温度'].values
interpolated_temps = temp_field[~np.isnan(temp_field)]
# 绘制直方图
ax.hist(measured_temps, bins=15, alpha=0.7, color='blue',
label=f'实测温度 (n={len(measured_temps)})', density=True)
ax.hist(interpolated_temps, bins=30, alpha=0.5, color='red',
label=f'插值温度 (n={len(interpolated_temps)})', density=True)
ax.set_xlabel('温度 (°C)', fontsize=10)
ax.set_ylabel('频率密度', fontsize=10)
ax.set_title(f'温度分布直方图\n{self.format_time_point(time_point)}', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
def run_analysis(self):
"""
运行完整的温度场分析
"""
print("开始边坡温度场分析...")
# 1. 生成每个时间点的温度场等值线云图
print("\n1. 生成每个时间点的温度场等值线云图")
for time_point in self.time_points:
self.create_time_temperature_contour(time_point)
# 2. 生成不同深度温度分析
print("\n2. 生成不同深度温度场分析")
self.create_depth_temperature_analysis()
# 3. 生成温度时间序列
print("\n3. 生成温度时间序列分析")
self.generate_temperature_time_series()
# 4. 生成分析报告
print("\n4. 生成分析报告")
self.generate_analysis_report()
print(f"\n分析完成！所有结果已保存到: {self.output_dir}")
def create_depth_temperature_analysis(self):
"""
创建不同深度的温度场分析 - 使用等值线
"""
print("正在生成不同深度的温度场分析...")
# 选择几个关键深度
depths = [self.z_grid[0],
self.z_grid[len(self.z_grid)//4],
self.z_grid[len(self.z_grid)//2],
self.z_grid[3*len(self.z_grid)//4],
self.z_grid[-1]]
# 选择几个关键时间点
num_times = min(3, len(self.time_points))
selected_times = [self.time_points[i * len(self.time_points) // num_times]
for i in range(num_times)]
# 为每个深度创建等值线图
for depth in depths:
fig, axes = plt.subplots(1, num_times, figsize=(6*num_times, 5))
if num_times == 1:
axes = [axes]
depth_index = np.argmin(np.abs(self.z_grid - depth))
for i, time_point in enumerate(selected_times):
temp_field = self.interpolate_temperature_3d(time_point)
if temp_field is None:
continue
# 在该深度的XY切片
slice_temp = temp_field[:, :, depth_index]
# 绘制等值线图（只绘制有效数据）
valid_mask = ~np.isnan(slice_temp)
if np.any(valid_mask):
contour = axes[i].contourf(self.x_grid, self.y_grid, slice_temp.T,
levels=20, cmap='jet')
# 添加等值线
contour_lines = axes[i].contour(self.x_grid, self.y_grid, slice_temp.T,
levels=10, colors='black', linewidths=0.5, alpha=0.7)
axes[i].clabel(contour_lines, inline=True, fontsize=6, fmt='%.1f°C')
# 添加监测点
time_data = self.merged_df[self.merged_df['时间'] == time_point]
scatter = axes[i].scatter(time_data['X坐标'], time_data['Y坐标'],
c=time_data['温度'], cmap='jet', s=30,
edgecolors='white', linewidth=0.5)
# 设置坐标轴范围，精确匹配模型范围
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
bounds = self.slope_mesh.bounds
axes[i].set_xlim(bounds[0][0], bounds[1][0])
axes[i].set_ylim(bounds[0][1], bounds[1][1])
axes[i].set_xlabel('X坐标 (mm)', fontsize=9)
axes[i].set_ylabel('Y坐标 (mm)', fontsize=9)
# 使用新的时间格式
formatted_time = self.format_time_point(time_point)
axes[i].set_title(f'{formatted_time}\n深度: {depth:.1f}mm', fontsize=10)
axes[i].grid(True, alpha=0.3)
plt.tight_layout()
# 添加共享颜色条
if np.any(valid_mask):
cbar = fig.colorbar(contour, ax=axes, shrink=0.8, orientation='horizontal',
pad=0.05, label='温度 (°C)')
# 保存图片
filename = f"深度_{depth:.1f}mm_温度等值线.png"
filepath = os.path.join(self.output_dir, filename)
plt.savefig(filepath, dpi=300, bbox_inches='tight')
plt.close()
print(f"已保存: {filename}")
def generate_temperature_time_series(self):
"""
生成温度时间序列分析 - 使用月-日 小时格式
"""
print("正在生成温度时间序列分析...")
# 为每个测点创建温度时间序列
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# 1. 所有测点的温度时间序列
for point_id in self.monitoring_points['测点编号'].unique():
point_data = self.merged_df[self.merged_df['测点编号'] == point_id]
point_data = point_data.sort_values('时间')
# 格式化时间为月-日 小时格式
formatted_times = point_data['时间'].apply(lambda x: self.format_time_point(x))
ax1.plot(formatted_times, point_data['温度'],
label=f'测点{point_id}', linewidth=1, alpha=0.7)
ax1.set_xlabel('时间 (月-日 小时)', fontsize=12)
ax1.set_ylabel('温度 (°C)', fontsize=12)
ax1.set_title('各测点温度时间序列', fontsize=14)
ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
ax1.grid(True, alpha=0.3)
# 如果时间点太多，旋转x轴标签避免重叠
if len(self.time_points) > 10:
plt.setp(ax1.get_xticklabels(), rotation=45, ha='right')
# 2. 平均温度时间序列
if '时间' in self.merged_df.columns:
# 按时间点计算平均温度
time_avg = self.merged_df.groupby('时间')['温度'].mean()
# 格式化时间为月-日 小时格式
formatted_times = [self.format_time_point(t) for t in time_avg.index]
ax2.plot(formatted_times, time_avg.values, 'ro-', linewidth=2, markersize=4)
ax2.set_xlabel('时间 (月-日 小时)', fontsize=12)
ax2.set_ylabel('平均温度 (°C)', fontsize=12)
ax2.set_title('时间点平均温度变化', fontsize=14)
ax2.grid(True, alpha=0.3)
# 如果时间点太多，旋转x轴标签避免重叠
if len(time_avg) > 10:
plt.setp(ax2.get_xticklabels(), rotation=45, ha='right')
plt.tight_layout()
# 保存图片
filename = "温度时间序列分析.png"
filepath = os.path.join(self.output_dir, filename)
plt.savefig(filepath, dpi=300, bbox_inches='tight')
plt.close()
print(f"已保存: {filename}")
def generate_analysis_report(self):
"""
生成分析报告
"""
# 基本统计信息
temp_stats = {
'最高温度': self.merged_df['温度'].max(),
'最低温度': self.merged_df['温度'].min(),
'平均温度': self.merged_df['温度'].mean(),
'温度标准差': self.merged_df['温度'].std()
}
coord_stats = {
'X坐标范围': f"{self.coord_df['X坐标'].min():.1f} - {self.coord_df['X坐标'].max():.1f} mm",
'Y坐标范围': f"{self.coord_df['Y坐标'].min():.1f} - {self.coord_df['Y坐标'].max():.1f} mm",
'高程范围': f"{self.coord_df['高程'].min():.1f} - {self.coord_df['高程'].max():.1f} mm"
}
if hasattr(self, 'slope_mesh') and self.slope_mesh is not None:
mesh_info = f"""
模型信息:
- 顶点数: {len(self.slope_mesh.vertices)}
- 面片数: {len(self.slope_mesh.faces)}
- 模型边界: {self.slope_mesh.bounds}
- 模型是否水密: {self.is_watertight}
"""
else:
mesh_info = "- 模型: 基于测点创建的简化模型"
report_content = f"""
边坡温度场分析报告
生成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}
数据概况:
- 数据时间范围: {self.format_time_point(self.time_points[0]) if self.time_points else '未知'} 至 {self.format_time_point(self.time_points[-1]) if self.time_points else '未知'}
- 总时间点: {len(self.time_points)}
- 监测点数量: {len(self.monitoring_points)}
- 温度记录数: {len(self.merged_df)}
温度统计:
- 最高温度: {temp_stats['最高温度']:.2f} °C
- 最低温度: {temp_stats['最低温度']:.2f} °C
- 平均温度: {temp_stats['平均温度']:.2f} °C
- 温度标准差: {temp_stats['温度标准差']:.2f} °C
坐标范围:
- X坐标: {coord_stats['X坐标范围']}
- Y坐标: {coord_stats['Y坐标范围']}
- 高程: {coord_stats['高程范围']}
切片位置设置:
- XY平面Z坐标: {self.xy_slice_z if self.xy_slice_z is not None else '自动(平均高度)'} mm
- XZ平面Y坐标: {self.xz_slice_y if self.xz_slice_y is not None else '自动(中心位置)'} mm
- YZ平面X坐标: {self.yz_slice_x if self.yz_slice_x is not None else '自动(中心位置)'} mm
{mesh_info}
生成文件说明:
1. 温度场等值线_YYYYMMDD_HHMMSS.png - 每个时间点的多视图温度场等值线图
2. 深度_XX.Xmm_温度等值线.png - 不同深度的温度等值线
3. 温度时间序列分析.png - 测点温度变化趋势
文件保存位置: {self.output_dir}
"""
report_path = os.path.join(self.output_dir, "分析报告.txt")
with open(report_path, 'w', encoding='utf-8') as f:
f.write(report_content)
print("已保存: 分析报告.txt")
# ========================= 示例使用 =========================
if __name__ == "__main__":
# 使用示例 - 只需要提供文件名，完整路径会自动构建
# 可以自定义切片位置（单位：mm）
analyzer = SlopeTemperatureAnalysis(
temp_data_file='温度数据.xlsx', # 文件位于 D:\桌面\温度数据\三维温度图\温度数据.xlsx
coord_data_file='坐标数据.xlsx', # 文件位于 D:\桌面\温度数据\三维温度图\坐标数据.xlsx
stl_file='complex_model.stl', # 文件位于 D:\桌面\温度数据\三维温度图\complex_model.stl
xy_slice_z=300, # XY平面切片的Z坐标(mm)
xz_slice_y=250, # XZ平面切片的Y坐标(mm)
yz_slice_x=300 # YZ平面切片的X坐标(mm)
)
# 运行分析
analyzer.run_analysis() 分析这个代码，解释每一步的含义，可以在哪里改动
最新发布
扫一扫
6482
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
