标识第一个图层中穿过所有3d numpy Z图层中高亮显示的属性区域的格网像元

6qqygrtg  于 2023-08-05  发布在  其他
关注(0)|答案(1)|浏览(87)

我有一个3D NumPy阵列,代表地下储层的孔隙度值,由40层组成,孔隙度分布各不相同。每个层被划分成大小为50 x50的网格,从而产生形状为(40,50,50)的孔隙率阵列,称为“porosity_array”。我想确定第一层中可能穿过所有层中特定孔隙度区域的网格单元。
为此,我将孔隙度值分为三类:贫穷,美丽,善良。分类范围定义如下:

poor_class = (np.min(porosity_array), 0.1)
    fair_class = (0.1, 0.14)
    good_class = (0.14, np.max(porosity_array))

字符串
或者我们可以将范围分为5类:非常差(np.min(porosity_array),0.1)、差(0.1,0.13)、一般(0.13,0.16)、好(0.16,0.19)和非常好(0.19,np.max(porosity_array))。
这种选择显然假设我们在每个类别中具有足够的数据点,使得当提出井位置时,每个孔隙度类别具有代表性的井位置。
使用具有“射流”色图的热图,红色表示高孔隙率,蓝色表示差孔隙率,黄色表示中等孔隙率。这些颜色基于孔隙率分布针对每个层而变化。
我想找到第一层中的网格单元,这些网格单元落在高孔隙度范围(红色区域)内,其中垂直红线将主要穿过所有储层的红色区域。同样,我的目标是识别第一层中的网格单元,用于低孔隙率范围(蓝色区域)和中等孔隙率范围(黄色区域),这些网格单元与所有层中的相应颜色对齐。显然,一个人不能穿透一个井在所有的层和交叉只有红色。所以越接近满足这个条件越好。使用较小的西格玛值,我们将加宽孔隙率的范围,但其代价是增加的随机性和异质性以及不同颜色上相同颜色的较小批次,因此使其穿过所有层上相同颜色(区域)的机会较小。我认为获胜的解决方案将是一个捕捉这些类别的孔隙率之间的过渡区,试图增加条件生存的机会。以下是我设想的识别这些网格单元的过程:
1.使用NumPy的where函数确定第一层中落在高孔隙率范围(红色区域)内的网格单元的指数。
1.迭代其余层(不包括第一层),并检查相同索引处的相应网格单元是否落入高孔隙率范围内。
1.计算满足条件的层数。
1.对低孔隙率范围(蓝色区域)和中等孔隙率范围(黄色区域)重复相同的过程,相应地调整条件。通过对满足条件的层的数量进行计数,我们可以识别第一层中更有可能穿过所有储层的对应颜色区域的网格单元。所有层中的垂直井必须相同,并且无论层号如何,垂直井必须保持在相同的网格单元中。建议的威尔斯也必须尽可能分散,而不是仅集中在储层的一个区域。
平滑数组是为了使所提供的数据尽可能真实地反映我自己的数据。
将非常感谢任何有效地完成这项任务的建议或替代方法。我下面的代码给出了所有三个区域的0层!

import numpy as np
    import matplotlib.pyplot as plt
    from scipy.ndimage import gaussian_filter

# Step 1: Create an empty 3D NumPy array
    porosity_array = np.zeros((40, 50, 50))

# Step 2: Generate random values within the desired range
    np.random.seed(123)
    porosity_array = np.random.uniform(low=0, high=0.25, size=(40, 50, 50))

# Step 3: Define the ranges for the three modes
   poor_class = (np.min(porosity_array), 0.1)
   fair_class = (0.1, 0.14)
   good_class = (0.14, np.max(porosity_array))

# Step 4: Assign values within each mode's range
   porosity_array[(porosity_array >= poor_class[0]) & (porosity_array < poor_class[1])] = 0.05  # Poor mode
   porosity_array[(porosity_array >= fair_class[0]) & (porosity_array < fair_class[1])] = 0.12  # Fair mode
   porosity_array[(porosity_array >= good_class[0]) & (porosity_array <= good_class[1])] = 0.2  # Good mode

# Apply Gaussian smoothing to make regions more continuous
   sigma = 1.0  # Lower sigma means wider ranges of porosity, although at a cost of having more randomness and heterogeneity and at a cost of smaller batches of same color across different colors.  Adjust the sigma value for desired smoothing effect
   porosity_array_smoothed = gaussian_filter(porosity_array, sigma=sigma)

# Plotting the heatmap
   plt.imshow(porosity_array_smoothed[0], cmap='jet', origin='lower')
   plt.colorbar(label='Porosity')
   plt.title('Porosity Heatmap (Smoothed)')
   plt.show()


   porosity_array = porosity_array_smoothed.copy()

# Define the porosity ranges

    poor_class = (np.min(porosity_array), 0.1)
    fair_class = (0.1, 0.14)
    good_class = (0.14, np.max(porosity_array))
    
    # Find the indices of grid cells in the first layer that fall within the high porosity range
    red_indices = np.where((porosity_array[0] >= good_class[0]) & (porosity_array[0] < good_class[1]))
    
    # Count the number of layers in which the condition is satisfied
    red_count = 0
    for layer in porosity_array[1:]:
        if np.all(layer[red_indices] >= good_class[0]) and np.all(layer[red_indices] < good_class[1]):
            red_count += 1
    
    # Find the indices of grid cells in the first layer that fall within the low porosity range
    blue_indices = np.where((porosity_array[0] >= poor_class[0]) & (porosity_array[0] < poor_class[1]))
    
    # Count the number of layers in which the condition is satisfied
    blue_count = 0
    for layer in porosity_array[1:]:
        if np.all(layer[blue_indices] >= poor_class[0]) and np.all(layer[blue_indices] < poor_class[1]):
            blue_count += 1
    
    # Find the indices of grid cells in the first layer that fall within the fair porosity range
    yellow_indices = np.where((porosity_array[0] >= fair_class[0]) & (porosity_array[0] < fair_class[1]))
    
    # Count the number of layers in which the condition is satisfied
    yellow_count = 0
    for layer in porosity_array[1:]:
        if np.all(layer[yellow_indices] >= fair_class[0]) and np.all(layer[yellow_indices] < fair_class[1]):
            yellow_count += 1
    
    # Print the results
    print("High Porosity (Red) - First Layer Indices:", red_indices)
    print("Number of Layers with High Porosity (Red):", red_count)
    print("Low Porosity (Blue) - First Layer Indices:", blue_indices)
    print("Number of Layers with Low Porosity (Blue):", blue_count)
    print("Fair Porosity (Yellow) - First Layer Indices:", yellow_indices)
    print("Number of Layers with Fair Porosity (Yellow):", yellow_count)


x1c 0d1x的数据

cngwdvgl

cngwdvgl1#

因此,我假设路径应该排序,使得最佳路径具有最多的红色层,红色间隙在可能的情况下由黄色填充,如果没有黄色,则为蓝色。
要做到这一点,我会数多少红色,黄色和蓝色层在50x50网格中的每个点。然后,我将按最大红色,然后黄色,然后蓝色对结果进行排序。

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt

plt.close("all")

porosity_array = np.zeros((40, 50, 50))

np.random.seed(123)
porosity_array = np.random.uniform(low=0, high=0.25, size=(40, 50, 50))

poor_class = (np.min(porosity_array), 0.1)
fair_class = (0.1, 0.14)
good_class = (0.14, np.max(porosity_array))

red = (porosity_array >= good_class[0]) & (porosity_array < good_class[1])
yellow = (porosity_array >= fair_class[0]) & (porosity_array < fair_class[1])
blue = (porosity_array >= poor_class[0]) & (porosity_array < poor_class[1])

porosity_array[blue] = 0.05  # Poor mode
porosity_array[yellow] = 0.12  # Fair mode
porosity_array[red] = 0.2  # Good mode

sigma = 1.0  # Adjust the sigma value for desired smoothing effect
porosity_array = gaussian_filter(porosity_array, sigma=sigma)

def get_porosity_counts(porosity_array, poor_class, fair_class, good_class):
    red = (porosity_array >= good_class[0]) & (porosity_array < good_class[1])
    yellow = (porosity_array >= fair_class[0]) & (porosity_array < fair_class[1])
    blue = (porosity_array >= poor_class[0]) & (porosity_array < poor_class[1])

    red_counts = np.sum(red, axis=0)
    # ignore locations already counted as red
    yellow[red] = False
    yellow_counts = np.sum(yellow, axis=0)
    # ignore locations already counted as red or yellow
    blue[red | yellow] = False
    blue_counts = np.sum(blue, axis=0)

    return red_counts, yellow_counts, blue_counts

def get_best_wells(red_counts, yellow_counts, blue_counts, n_best=None):
    # sort by largest red count, then yellow, then blue
    # reverse because the sort is from smallest to largest
    best_wells = np.lexsort((blue_counts.ravel(),
                             yellow_counts.ravel(),
                             red_counts.ravel()))[::-1]
    best_wells_indices = np.unravel_index(best_wells, red_counts.shape)
    best_wells_indices = [(i, j) for i, j in zip(best_wells_indices[0],
                                                 best_wells_indices[1])]

    if n_best:
        best_wells_indices = best_wells_indices[:n_best]

    return best_wells_indices

def get_metrics(index, red_counts, yellow_counts, blue_counts):
    red = red_counts[index]
    yellow = yellow_counts[index]
    blue = blue_counts[index]

    return {"red": red, "yellow": yellow, "blue": blue, "index":index}

poor_class = (np.min(porosity_array), 0.1)
fair_class = (0.1, 0.14)
good_class = (0.14, np.max(porosity_array))

red = (porosity_array >= good_class[0]) & (porosity_array < good_class[1])
yellow = (porosity_array >= fair_class[0]) & (porosity_array < fair_class[1])
blue = (porosity_array >= poor_class[0]) & (porosity_array < poor_class[1])

counts = get_porosity_counts(porosity_array, poor_class, fair_class, good_class)
indices = get_best_wells(*counts)
print(get_metrics(indices[0], *counts))  # {'red': 19, 'yellow': 18, 'blue': 3, 'index': (0, 0)}

fig, ax = plt.subplots(figsize=(10,10))
x = np.arange(porosity_array[0].shape[0])
y = np.arange(porosity_array[0].shape[1])
ax.pcolormesh(x, y, porosity_array[0], cmap="jet")
for i, index in enumerate(get_best_wells(*counts, 10), start=1):
    ax.text(*index, i, ha="center", va="center")
ax.set_aspect(1)
fig.show()

字符串


的数据
下面是我们改变sigma=2.0的结果:


相关问题