如何将需要插值的python代码矢量化为特定的数据点

2024-06-16 12:38:56 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个问题,我用一个叫做MCNP的计算机程序,从一个粒子通量计算正方形几何体中的能量沉积。正方形几何体被分解为一个网格网格,网格的长度、宽度和高度为50个立方网格。数据被放入一个文本文件中,显示每个网格在笛卡尔坐标系(x,y,z位置)的质心位置以及在该x,y,z坐标系下的能量沉积。然后使用Python脚本提取数据。我有一个脚本,可以让我在z平面上切一片,在这个平面上绘制能量沉积的热图,这个脚本可以工作,但是我认为它不是很有效,我正在寻找将过程矢量化的解决方案。在

该代码读取X、Y和Z坐标作为三个独立的1-D numpy阵列,并将该坐标处的能量沉积读取为1-D numpy阵列。为了说明这一点,我们假设我想在零的Z坐标处取一个切片,但是没有一个网格质心位于0的Z坐标,然后,我必须(并且确实)遍历Z坐标数组中的所有数据点,直到它找到一个大于0(数组索引I)的数据点,并且数组索引(I-1)小于0。然后,它需要使用Z空间中的这些阵列点以及切片位置(在本例中为0)和这些阵列指数处的能量沉积,并进行插值,以在切片的Z位置找到正确的能量沉积。因为X和Y阵列不受影响,现在我有了X,Y的坐标,可以绘制出特定X,Y位置的热图和切片位置的能量沉积。代码还需要确定切片位置是否已经在数据集中,在这种情况下不需要插值。我的代码可以工作,但是我不知道如何使用内置的scipy插值方案,而是编写了一个函数来执行插值。在这个场景中,必须使用for循环进行迭代,直到找到z位置位于切片位置的上方和下方的位置(本例中z=0)。我在这篇文章中附上了我的示例代码,并请求帮助更好地矢量化这个代码片段(如果它可以更好地矢量化),并希望在这个过程中学到一些东西。在

# - This transforms the read in data from a list to a numpy array
#   where Magnitude represents the energy deposition
XArray = np.array(XArray); YArray = np.array(YArray)
ZArray = np.array(ZArray); Magnitude = np.array(Magnitude)

#==============================================================
# - This section creates planar data for a 2-D plot

# Interpolation function for determining 2-D slice of 3-D data
def Interpolate(X1,X2,Y1,Y2,X3):
    Slope = (Y2-Y1)/(X2-X1)
    Y3 = (X3-X1)*Slope
    Y3 = Y3 + Y1
    return Y3

# This represents the location on the Z-axis where a slice is taken
Slice_Location = 0.0

XVal = []; YVal = []; ZVal = []
Tally = []; Error = []
counter = 1.0
length = len(XArray)-1
for numbers in range(length):
 #   - If data falls on the selected plane location then use existing data
    if ZArray[counter] == Slice_Location:       
        XVal.append(XArray[counter])
        YVal.append(YArray[counter])
        ZVal.append(ZArray[counter])
        Tally.append(float(Magnitude[counter]))

# - If existing data does not exist on selected plane then interpolate
    if ZArray[counter-1] < Slice_Location and ZArray[counter] > Slice_Location:
        XVal.append(XArray[counter])
        YVal.append(YArray[counter])
        ZVal.append(Slice_Location)
        Value = Interpolate(ZArray[counter-1],ZArray[counter],Magnitude[counter-1], \
                           Magnitude[counter],Slice_Location)
        Tally.append(float(Value))
    counter = counter + 1

XVal = np.array(XVal); YVal = np.array(YVal); ZVal = np.array(ZVal)
Tally = np.array(Tally);

Tags: 数据代码网格datanpcounter切片slice