通过由一串点定义的横截面在固定网格上产生三维高程-Java 学习之路

我试图在1米或0.5米间距的固定网格上创建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米，有一个双梯形截面 .

我试过，编写这个但是没有约会：

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

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

码：

设置计算域

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 import 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（）

1 回答

0
您可以通过最初将第三个维度设置为零来从头开始处理您的渠道配置文件 . 通过这种方式，您将能够沿曲线旋转和平移您的轮廓 . 例如：
```
#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）“像素大小”，换句话说，您的 Channels 配置文件以米为单位，但在DEM中，您正在使用矩阵中的（基本上）索引 . 需要 Build 一个简单的线性变换，以便确定“距离中心线”-20米的“像素数”

2）Translate（）和Rotate（）函数 . 这里有任何简单的矢量数学 . 请注意，如果您通过参考0,0,0来表达您的 Channels 配置文件，则旋转将是非常简单的表达式 . 有关详细信息，请参阅此处：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的一个"sticking out" .

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

希望这有所帮助 . 您是在处理CFD问题还是仅用于可视化？
回复于 2024-05-09T07:55:38+08:00

通过由一串点定义的横截面在固定网格上产生三维高程

设置计算域

1 回答

相关问题