|
|
|
|
@ -0,0 +1,207 @@
|
|
|
|
|
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)
|
|
|
|
|
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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
data_name = ['data.pos', 'data1.pos']
|
|
|
|
|
level_nums = 15
|
|
|
|
|
gWeight = 16
|
|
|
|
|
data = create()
|
|
|
|
|
fig, ax = plt.subplots()
|
|
|
|
|
zmin=1
|
|
|
|
|
zmax=0
|
|
|
|
|
for i in range(len(data_name)):
|
|
|
|
|
data.data_pre(data_name[i])
|
|
|
|
|
data.contours_pre(level_nums, gWeight)
|
|
|
|
|
|
|
|
|
|
"""画出密度等高线"""
|
|
|
|
|
f = data.f
|
|
|
|
|
x_min, x_max = data.x_range[0], data.x_range[1]
|
|
|
|
|
y_min, y_max = data.y_range[0], data.y_range[1]
|
|
|
|
|
levels = data.levels
|
|
|
|
|
|
|
|
|
|
xx, yy = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]
|
|
|
|
|
# 填充等高线图中间的区域
|
|
|
|
|
cmap = mpl.colormaps.get_cmap('jet')
|
|
|
|
|
colors = cmap(np.linspace(0, 1, len(data_name)))
|
|
|
|
|
level = [levels[0], levels[-1]]
|
|
|
|
|
ff = ax.contourf(xx, yy, f, levels=level, colors=colors[len(data_name)-i-1:len(data_name)-i], alpha=1, zorder=len(data_name)-i)
|
|
|
|
|
zmin = min(zmin, ff.zmin)
|
|
|
|
|
zmax = max(zmax, ff.zmax)
|
|
|
|
|
plt.xlabel('X')
|
|
|
|
|
plt.ylabel('Y')
|
|
|
|
|
plt.title('Density Distribution')
|
|
|
|
|
|
|
|
|
|
# 生成colorbar
|
|
|
|
|
sm = mpl.cm.ScalarMappable(cmap=cmap, norm=mpl.colors.Normalize(vmin=zmin, vmax=zmax))
|
|
|
|
|
sm.set_array([])
|
|
|
|
|
cbar = fig.colorbar(sm, ax=ax)
|
|
|
|
|
|
|
|
|
|
wells_name = 'well.txt'
|
|
|
|
|
wells = well_to_edge()
|
|
|
|
|
wells.wells_name_and_position(wells_name)
|
|
|
|
|
wells.welltoedge_distance(data.contours)
|
|
|
|
|
|
|
|
|
|
"""画井位信息"""
|
|
|
|
|
typee = wells.type
|
|
|
|
|
points = wells.position
|
|
|
|
|
namee = wells.name
|
|
|
|
|
min_distance = wells.min_distance
|
|
|
|
|
contours_p = wells.welltoedge_points
|
|
|
|
|
for i in range(len(points)):
|
|
|
|
|
if typee[i][0] == 0:
|
|
|
|
|
ax.scatter(points[i][0], points[i][1], marker='o', edgecolors='black', facecolors='none', s=100)
|
|
|
|
|
ax.scatter(points[i][0], points[i][1], marker='o', edgecolors='black', facecolors='none', s=50,
|
|
|
|
|
linewidths=1)
|
|
|
|
|
ax.scatter(points[i][0], points[i][1], marker='o', edgecolors='black', facecolors='none', s=30)
|
|
|
|
|
else:
|
|
|
|
|
ax.scatter(points[i][0], points[i][1], marker='o', edgecolors='black', facecolors='black', s=100, zorder=len(data_name)+1)
|
|
|
|
|
ax.annotate(namee[i][0], (points[i][0], points[i][1]), xytext=(5, 5), textcoords='offset points')
|
|
|
|
|
|
|
|
|
|
# 绘制油井边缘最短距离点的箭头并标出距离
|
|
|
|
|
ax.quiver(points[i][0], points[i][1], contours_p[i][0] - points[i][0], contours_p[i][1] - points[i][1],
|
|
|
|
|
angles='xy', scale=1.03,
|
|
|
|
|
scale_units='xy', width=0.002, zorder=len(data_name)+1) # 绘制箭头
|
|
|
|
|
ax.text(points[i][0] / 2 + contours_p[i][0] / 2 + min_distance[i] / 20,
|
|
|
|
|
points[i][1] / 2 + contours_p[i][1] / 2 - min_distance[i] / 20, f'{round(min_distance[i], 2)}',
|
|
|
|
|
fontdict={'size': '10', 'color': 'g'}) # 标出距离
|
|
|
|
|
ax.text(points[i][0] * 3 / 4 + contours_p[i][0] / 4 + min_distance[i] / 20,
|
|
|
|
|
points[i][1] * 3 / 4 + contours_p[i][1] / 4 - min_distance[i] / 20,
|
|
|
|
|
r'$\theta$=' f'{round(wells.angle[i], 2)}', fontdict={'size': '10', 'color': 'y'}) # 标出角度
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
冯冬卫