You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

135 lines
4.7 KiB
Python

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

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