STL中sort的使用 _Updated 150721

本文详细介绍 STL 中 sort 函数的使用方法,包括基本用法、如何指定自定义比较函数实现不同类型数据的排序,以及通过具体示例展示如何应用于结构体数组的排序。

 

STLsort函数用法简介

ACM题的时候,排序是一种经常要用到的操作。如果每次都自己写个冒泡之类的O(n^2)排序,不但程序容易超时,而且浪费宝贵的比赛时间,还很有可能写错。STL里面有个sort函数,可以直接对数组排序,复杂度为n*log2(n)

使用这个函数,需要包含头文件#include <algorithm>

这个函数可以传两个参数或三个参数。第一个参数是要排序的区间首地址,第二个参数是区间尾地址的下一地址。也就是说,排序的区间是[a,b)。简单来说,有一个数组int a[100],要对从a[0]a[99]的元素进行排序,只要写sort(a,a+100)就行了,默认的排序方式是升序。

拿我出的AC的策略这题来说,需要对数组t的第0len-1的元素排序,就写sort(t,t+len);

对向量 v排序也差不多,sort(v.begin(),v.end());

排序的数据类型不局限于整数,只要是定义了小于运算的类型都可以,比如字符串类string

如果是没有定义小于运算的数据类型,或者想改变排序的顺序,就要用到第三参数——比较函数。比较函数是一个自己定义的函数,返回值是bool型,它规定了什么样的关系才是小于。想把刚才的整数数组按降序排列,可以先定义一个比较函数cmp

bool cmp(int a,int b)  { return a>b; } 


排序的时候就写 sort(a,a+100,cmp);

假设自己定义了一个结构体 node

struct node{ 
int a; 
int b; 
double c; 
} 


 

有一个 node类型的数组node arr[100],想对它进行排序:先按a值升序排列,如果a值相同,再按b值降序排列,如果b还相同,就按c降序排列。就可以写这样一个比较函数:

以下是代码片段:

bool cmp(node x,node y) { 
if(x.a!=y.a) return x.a>y.b;
if(x.b!=y.b) return x.b>y.b; 
return x.c>y.c; 
} 


 

排序时写 sort(arr,a+100,cmp);

最后看一个完整的实例,初赛时的一道题目文件名排序

以下是代码片段:

#include<iostream> 
#include<algorithm> 
#include<string> 
using namespace std; 
// 定义一个结构体来表示文件, a 代表文件名, b 代表文件类型(要么 "File" 要么 "Dir" )
struct node{ 
string a,b; 
}; 

//ASCII 码中,所有大写字母排在所有小写字母前面, 'A'<'Z'<'a'<'z' 
// 而这题要求忽略大小写,所以不能直接用字符串的比较。自定义了一个 lt 函数,就是 less than 的意思
// 先把两个字符串全部转化为小写,再比较大小(字典序)

bool lt(string x,string y) 
{ 
int i; 
for(i=0;i<x.length();i++) 
if(x[i]>='A'&&x[i]<='Z') 
x[i]='a'+(x[i]-'A'); 
for(i=0;i<y.length();i++) 
if(y[i]>='A'&&y[i]<='Z') 
y[i]='a'+(y[i]-'A'); 
return x<y; 
} 

// 自定义的比较函数,先按 b 值升序排列(也就是 "Dir" 排在 "File" 前面)
// 如果 b 值相同,再按 a 升序排列,用的是刚才定义的 lt 函数

bool comp(node x,node y) 
{ 
if(x.b!=y.b)return x.b<y.b; 
return lt(x.a,y.a); 
} 

int main() 
{ 
node arr[10001]; 
int size=0; 
while(cin>>arr[size].a>>arr[size].b) 
size++; 
sort(arr,arr+size,comp); 
for(int i=0;i<size;i++) 
cout<<arr[i].a<<" "<<arr[i].b<<endl; 
return 0; 
} 


<Updated 20150721>

感谢Zoecur指出上述代码中的两处错误,已改正,为两年前的不注意道歉。

以及,之前在竞赛中遇到过sort和qsort的混淆问题,弄清之后也贴于此处。 

sort和qsort所需要的cmp函数不同,并且调用的时候的姿势也不同——


【sort_cmp.cpp】

#include <queue>
#include <cstdio>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

#define MAXN 100

struct In 
{    
	int data;   
	char str[100];   
}s[MAXN];   
  
// 按照结构体中字符串str字典序排序   
inline int cmp (const In a, const In b){return strcmp( a.str , b.str );}

struct Point  
{   
	int x;   
	int y;   
}p[MAXN];   
  
// 按照 x 从小到大排序,当 x 相等时按照 y 从大到小排序   
int cmp2( const Point a , const Point b )
{   
	if(a.x != b.x) return a.x - b.x;   
	else return b.y - a.y;   
}   

int main()
{
	sort(s,s+MAXN,cmp);  
	sort(p,p+MAXN,cmp2); 
	return 0;	
} 


【qsort_cmp.cpp】
#include <queue>
#include <cstdio>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

#define MAXN 100

struct In 
{    
	int data;   
	char str[100];   
}s[MAXN];   
  
// 按照结构体中字符串str字典序排序   
int cmp ( const void *a , const void *b )   
{   
	return strcmp( (*(In *)a).str , (*(In *)b).str );   
}   

struct Point  
{   
	int x;   
	int y;   
}p[MAXN];   
  
// 按照 x 从小到大排序,当 x 相等时按照 y 从大到小排序   
int cmp2( const void *a , const void *b )
{   
	struct Point *c = (Point *)a;   
	struct Point *d = (Point *)b;   
	if(c->x != d->x) return c->x - d->x;   
	else return d->y - c->y;   
}   

int main()
{
	qsort(s,MAXN,sizeof(s[0]),cmp);  
	qsort(p,MAXN,sizeof(p[0]),cmp2); 
	return 0;	
} 

关于较为常用的凸包排序cmp函数,今年校赛也有设计到,将此题也贴于此处方便查阅

校赛Solution传送门

Problem_[G]raphic

timelimit 2 seconds / memory limit 256 megabytes

 

Description

Curabis always do excellent job on computation geometry problems. Indeed, he has acquired the first solution award about the computation geometry problems in the ACM-ICPC Asia Regional several times. However, JKi, his stupid teammate, do it much worse. 

In order to help poor JKi, Curabis decide provide him with some practice problem set about computation geometry.

One of those problems is this:

Give a set of N2d-points on the plane, is there any three 2d-points of the 2d-point set areable to constructive an triangle, and all 2d-points of this set are either inside of the triangle or on the edge of the triangle?

Would you please help JKi finishing this program?

 

Input

The first line of the INPUT contain a positive integer T(T<=20) which denotes the number of test cases.

The first line of each case contains an integer N.

The next N(3<=N<=1000) lines of each case contain the 2d-points set. The i-th line contains two real numbers Xi, Yi denote the position of the 2d-point.

It is guaranteed that -1000.00<= Xi, Yi <=1000.00, and all real numbers in the input keep two decimal places exactly.

It is guaranteed that all 2d-points in any case can not stand on one straight line over all.

 

Output

In each line print an upper-case letter — ‘Y’ for the triangle is existed, otherwise ‘N’ instead.

 

Sample Input:

Sample Output:

3

4

0.00 0.00

0.00 2.00

2.00 0.00

0.50 0.50

4

0.00 0.00

0.00 2.00

2.00 0.00

2.00 2.00

5

0.00 0.00

0.00 0.00

1.00 0.00

2.00 0.00

1.00 1.00

Y

N

Y

 

梗: 其实我想说这个是不是博熊的ID啊…… <Update>啊啦被我猜对了~~

题意:有好多个点,给所有点的坐标,问是否存在三个点,他们构成的三角形可以包住所有的点(在线上和在点上也算包住)

思路:啊,包住,凸包嘛~ 题目等效于,这些点的凸包定点数是否不大于3,是的话输出Y,反之输出N

(可以作为凸包模板哟)


#include <cmath>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const double EPS=1e-10;

double add(double a,double b)
{
    if( fabs(a+b) < EPS*( abs(a)+abs(b)) )return 0;
    return a+b;
}

struct P
{
    double x;
    double y;
    P(){}
    P(double x,double y):x(x),y(y){}
    P operator+ (P p)
    {return P(add(x,p.x),add(y,p.y));}
    P operator- (P p)
    {return P(add(x,-p.x),add(y,-p.y));}
    P operator* (double d)
    {return P(x*d,y*d);}
    double dot(P p)//neiji
    {return add(x*p.x,y*p.y);}
    double det(P p)//waiji
    {return add(x*p.y,-y*p.x);}
};

bool cmp_x(const P& p,const P& q)
{
    if(p.x != q.x) return p.x<q.x;
    return p.y<q.y;
}

vector<P> convex_hull(P* ps,int n)
{
    sort(ps,ps+n,cmp_x);
    //for(int i=0;i<n;i++){cout<<ps[i].x<<":"<<ps[i].y<<endl;}
    int k=0;
    vector<P> qs(n<<1);
    for(int i=0;i<n;i++)
    {
        while( k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0 )k--;
        qs[k++]=ps[i];
    }
    for(int i=n-2,t=k; i>=0; i--)
    {
        while( k>t && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0 )k--;
        qs[k++]=ps[i];
    }
    qs.resize(k-1);
    return qs;
}

int N;
P ps[1024];

void solve()
{
    vector<P> qs = convex_hull(ps,N);
    //cout<<qs.size()<<endl;
    if(qs.size()<=3)printf("Y\n");
    else printf("N\n");
}

int main()
{
    freopen("G.in","r",stdin);
    int T; scanf("%d",&T);
    for(int _T=0;_T<T;_T++)
    {
        double tx,ty;
        scanf("%d",&N);
        memset(ps,0,sizeof ps);
        for(int i=0;i<N;i++)
        {
            cin>>tx>>ty;
            //cout<<tx<<":"<<ty<<endl;
            ps[i].x=tx;
            ps[i].y=ty;
        }
        solve();
    }
    return 0;
}




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() 分析这个代码,解释每一步的含义,可以在哪里改动
最新发布
10-16
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

糖果天王

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值