СОЗДАТЬ 3D-отметки над фиксированной сеткой из поперечного сечения, заданного строкой точек - PullRequest
Новая
       1

СОЗДАТЬ 3D-отметки над фиксированной сеткой из поперечного сечения, заданного строкой точек

0 голосов
/ 30 марта 2012

Я пытаюсь создать трехмерную поверхность на фиксированной сетке, скажем, с интервалом 1 м или 0,5 м, где поверхность - это канал, определяемый поперечным сечением нескольких точек. В идеале это должно быть любое количество баллов. Например, сечение, такое как:

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, тем самым предоставляя значения XYZ 3D ландшафта

В этом примере предполагается создать канал длиной 1500 м и шириной 70 м

КОД:

Настройка вычислительной области 

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 обеспечивает наклонное основание, а затем определенное поперечное сечение создает Z в диапазоне Y

В идеале это также может быть применено вдоль ломаной линии, а не только в направлении x. Таким образом, канал может быть получен вдоль кривой или S-изгибом, например

Надеюсь, у кого-то есть какие-нибудь умные идеи, как решить эту проблему ... Спасибо

Было упомянуто, что Сципи может помочь здесь ... Я попытаюсь разобраться в этом, чтобы создать функцию для интерполяции между точками:

из 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)

import matplotlib.pyplot as plt

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

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

plt.show ()

1 Ответ

0 голосов
/ 30 марта 2012

Вы можете обработать свой профиль канала в 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, которая является "торчащей" из вашей матрицы высот.

4) Ориентация ЦМР на реальный мир. В приведенном выше примере мы начинаем с 0.0 и будем двигаться на юг, а затем на юго-восток, потому что индексы в матрице увеличиваются сверху вниз

Надеюсь, это немного поможет. Вы имеете дело с проблемой CFD или это только для визуализации?

...