0

我正在尝试在间距为 1m 或 0.5m 的固定网格上创建 3D 表面，其中表面是由多个点的横截面定义的通道。理想情况下，它应该是任意数量的点。例如横截面，例如：

PTS = [[0.0,10.0],[3.0,9.0],[30.0,8.5],[33.0,8.0],[35.0,7.8],[37.0,8.0],[40.0,8.5],[67.0,9.0] ,[70.0,10.0]] 这里的河道宽 70 m，截面为双梯形。

我已经尝试过对此进行编码，但未能成功：

我想读取这些点，然后根据间距进行插值以提供计算出的高程（Z 值）。这应该填充 X 和 Y 的域，从而提供 3D 地形的 XYZ 值

这个例子应该创建一个 70 m 宽的 1500 m 长的通道

代码：

设置计算域

length = 50.0  # change back to 3500
width  = 30.0  #changed this
dx = dy = 1.0          # Resolution: of grid on both axes
h=2.21 # Constant depth
slope_X=1/100
def topography(x,y):
z = -x*slope_X
PTS = [[0.0,10.0],[3.0,9.0],[30.0,8.5],[33.0,8.0],[35.0,7.8],[37.0,8.0],[40.0,8.5],[67.0,9.0],[70.0,10.0]]


N = len(x)
for i in range(N):
    # Construct Cross section from LIST of PTS
    for j in range(len(PTS)):
        if j == 0:
            pass
        else:
            m = (PTS[j][1]-PTS[j-1][1])/(PTS[j][0]-PTS[j-1][0])
            b = PTS[j-1][1]-(PTS[j][1]-PTS[j-1][1])/(PTS[j][0]-PTS[j-1][0])*PTS[j-1][0]
            z[i]= m *y[i]+b

    if x[i]==10:
        print 'Z =', z[i]

return z

当代码遍历 X 时，基本 Z 提供了一个倾斜的床，然后定义的横截面在 Y 的范围内创建了 Z

理想情况下，这也可以沿多段线应用，而不是仅在 x 方向上应用。这样，通道可以沿曲线或 S 形弯曲产生

我希望有人对如何解决这个问题有一些聪明的想法......谢谢

有人提到 scipy 可能会在这里提供帮助....我将尝试理解这一点，以创建一个在点之间进行插值的函数：

从 scipy.interpolate 导入 interp1d

x = np.linspace(0, 10, 10)

y = np.exp(-x/3.0)

f = interp1d(x, y)

f2 = interp1d(x, y, kind='cubic')

xnew = np.linspace(0, 10, 40)

将 matplotlib.pyplot 导入为 plt

plt.plot(x,y,'o',xnew,f(xnew),'-',xnew, f2(xnew),'--')

plt.legend(['data', 'linear', 'cubic'], loc='best')

plt.show()

4

1 回答 1

0

您可以通过最初将第三维设置为零来从一开始就以 3D 方式处理您的频道配置文件。通过这种方式，您将能够沿曲线旋转和平移您的个人资料。例如：

#The DEM
DEM = numpy.array((height,width)); #...because a row corresponds to the y dimension
#Channel center line
cCurve = [[0.0,0.0],[0.0,1.0],[1.0,2.0]] #A channel going north and then turning North-East
#Channel profile. It is better if you express this in relative coordinates to the center line. In this case, the points left to the center line would have negative X values and the points to the right would have positive X values. 
PTS = [[0.0,0.0,10.0],[3.0,0.0,9.0],[30.0,0.0,8.5],[33.0,0.0,8.0],[35.0,0.0,7.8],[37.0,0.0,8.0],[40.0,0.0,8.5],[67.0,0.0,9.0],[70.0,0.0,10.0]];
#
for (aCenterLinePoint in range(0,len(cCurve)-1)):
     #Translate the channel profile to the current location of the center line
     translatedPTS = Translate(PTS,cCurve[aCenterLinePoint]);
     #Rotate the channel profile, along the Z axis to an angle that is defined by the current center line point and the next center line point
     rotatedTranslatedPTS = Rotate(rotatedPTS,getAngle(cCurve[aCenterLinePoint],cCurve[aCenterLinePoint+1]))
     # "Carve" the DEM with the Channel profile. You can apply interpolation here if you like
     for (aChannelPoint in rotatedTranslatedPTS):
         DEM[aChannelPoint[1], aChannelPoint[0]] = aChannelPoint[2]; #Please note the reversal of the X,Y coordinates to account for the classical array indexing!

我希望上面的片段传达了这个想法:-)。缺少的东西，您必须根据您的问题进行计算：

1）“像素大小”，换句话说，您的通道配置文件以米表示，但在 DEM 中，您正在使用（基本上）矩阵中的索引。需要建立一个简单的线性变换，这样就可以确定“多少像素”是“离中心线-20米”

2) Translate() 和 Rotate() 函数。任何简单的向量数学都可以在这里完成。请注意，如果您参考 0,0,0 来表达您的频道配置文件，则旋转将是非常简单的表达式。有关更多信息，请参见此处：http ://en.wikipedia.org/wiki/Rotation_matrix

3) getAngle() 函数是一个简单的 atan2(vectorA, vectorB)。（例如：http ://docs.scipy.org/doc/numpy/reference/generated/numpy.arctan2.html ）。请注意，在这种情况下，您将围绕 Z 轴旋转，这是 DEM 的“突出”轴。

4) DEM 对现实世界的定位。在上面的示例中，我们从 0.0 开始，我们将向南移动，然后向东南移动，因为矩阵中的索引自上而下增加

希望这有所帮助。您是在处理 CFD 问题还是只是为了可视化？

于 2012-03-30T11:22:07.077 回答