import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from scipy.stats import gaussian_kde from scipy.spatial.distance import cdist import copy class create(): def __init__(self): """ data is the x and y coordinate of data.pos contours is the edge points f is the kernel density value levels is the levels for contours x_range is the range of x label y_range is the range of y label """ self.data = [] self.contours = [] self.f = [] self.levels = [] self.x_range = [] self.y_range = [] def data_pre(self, data_name): with open(data_name, 'r') as f: lines = f.readlines() L_en=len(lines) lines= lines[1:L_en] data = [] for line in lines: x, y,z,t = line.strip().split("\t") data.append(list(map(float, [x,y]))) data = np.array(data) self.data = data self.x_range = [min(data[:, 0]), max(data[:,0])] self.y_range = [min(data[:, 1]), max(data[:,1])] def singlefile(self): pass def contours_pre(self, level_nums, gWeight): x = self.data[:, 0] y = self.data[:, 1] # 使用scipy库中的gaussian_kde函数计算密度估计 kde = gaussian_kde(self.data.T) # 生成网格点坐标 xx, yy = np.mgrid[x.min():x.max():200j, y.min():y.max():200j] positions = np.vstack([xx.ravel(), yy.ravel()]) # 计算网格点上的密度估计值 f = np.reshape(kde(positions).T, xx.shape) # 绘制等高线图 levels = [] for i in range(level_nums): levels.append(f.max() - (level_nums - i - 1) * (f.max() - f.min()) / gWeight) self.levels = levels contours = plt.contour(xx, yy, f, levels=[levels[0], levels[level_nums - 1]], cmap='coolwarm', alpha=0) plt.close() self_contours = [] for i in range(len(contours.allsegs[0])): self_contours += contours.allsegs[0][i].tolist() for i in range(len(self_contours)): if self_contours[i] not in self.contours: self.contours += [self_contours[i]] self.f = f class well_to_edge(): def __init__(self): """ name is used to store the well names type is the types of the wells position is the coordinates of wells min_distance is the minimum distances between wells and edge welltoedge_points is the points responding to the min_distance angle is the angles between the shortest distance direction vector from the well to the edge and the positive direction of the y-axis during clockwise rotation; wells_num: the number of wells """ self.name = [] self.type = [] self.position = [] self.min_distance = [] self.welltoedge_points = [] self.angle = [] self.wells_num = 0 def wells_name_and_position(self, wells_name): # 读取井位信息 with open(wells_name, 'r') as f_j: j_ing = f_j.readlines() points = [] typee = [] namee = [] for line in j_ing: x, y, type, name = line.strip().split("\t") points.append(list(map(float, [x, y]))) typee.append(list(map(int, [type]))) namee.append(name.split('\n')) self.position = points self.name = namee self.type = typee self.wells_num = len(points) def welltoedge_distance(self, contour_points): min_distance = [] contours_p=[] angles=[] wells_num = self.wells_num points = self.position # 定义待计算距离的点 for i in range(wells_num): point = np.array([points[i][0], points[i][1]]) # 计算点到等高线上所有点之间的距离 distances = cdist(point.reshape(1, -1), contour_points) # 取距离的最小值 min_distance.append(distances.min()) # 记录最短距离对应的点 for ii in range(distances.size): if distances[0][ii] == min_distance[-1]: contours_p.append(contour_points[ii]) # 计算角度(y轴为正北方向) direct = contours_p[-1] - point if direct[0] < 0: angles.append(360-np.arccos(np.dot(direct, np.array(([0, 1]))) / np.linalg.norm( direct)) / np.pi * 180) else: angles.append(np.arccos(np.dot(direct, np.array(([0, 1]))) / np.linalg.norm( direct)) / np.pi * 180) self.min_distance = min_distance self.welltoedge_points = contours_p self.angle = angles