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 がここで役立つ可能性があると述べられていました.... ポイント間を補間する関数を作成するために、これを理解しようとします。

from 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()

4

1 に答える 1

0

最初に 3 番目の次元をゼロに設定することで、チャネル プロファイルを最初から 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 )。この場合、DEM から「突き出ている」Z 軸を中心に回転することに注意してください。

4) 現実世界に対する DEM の向き。上記の例では、0.0 から開始し、マトリックス内のインデックスが上から下に増加するため、南に移動し、次に南東に移動します。

これが少し役立つことを願っています。CFDの問題を扱っていますか、それとも単に視覚化するためですか?

于 2012-03-30T11:22:07.077 に答える