|
|
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 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 |