commit fab74af61c2eca3d68d362105334947b31fd0eba
Author: 冯冬卫 <FDW@zhejianglab.com>
Date:   Wed Mar 29 17:04:49 2023 +0800

    代码提交CREATE.py

diff --git a/CREATE.py b/CREATE.py
new file mode 100644
index 0000000..e838820
--- /dev/null
+++ b/CREATE.py
@@ -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()
+
+
+
+
+