Txgy.WFLEA/CREAT.py

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