CReate 3D elevations over a fixed Grid from a Cross section defined by a string of Points

Question

I am trying to create a 3D surface on a fixed grid of say 1m or 0.5m spacing, where the surface is a channel defined by a cross section of a number of points. Ideally it should be any number of points. For example a cross section such as:

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]] Where the channel here is 70 m wide and has a double trapezoidal section.

I have tried, to code this but failed to date:

I want to read the points, and then interpolate based on the spacing to provide a calculated elevation (Z value). This should populate the domain of X & Y, thereby providing XYZ values of the 3D terrain

This example is supposed to create a 1500m long channel that is 70 m wide

CODE:

Setup computational domain

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

As the code traverses over X, the basic Z provides a sloping bed, and then the cross section defined, creates the Z's over a range of Y's

Ideally this could also be applied along a polyline, rather than apply it only in the x -direction. That way the channel could be produced along a curve, or a S- bend for example

I hope some one has some clever ideas on how to resolve this... Thank you

It was mentioned that scipy might be able to help here.... I will try to make sense of this, to create a function to interpolate between the points:

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)

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

A_A A_A · Accepted Answer · 2012-03-30T11:22:07

You can treat your channel profile in 3D right from the beginning by setting the third dimension to zero initially. In this way you will be able to rotate and translate your profile along a curve. For example:

#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!

I hope that the above snippet conveys the idea :-) . Things that are missing and you would have to calculate depending on your problem are:

1) The "pixel size", in other words, your channel profile is expressed in meters but in the DEM you are working with (essentially) indexes in a matrix. A simple linear transformation needs to be established so that you can determine "how many pixels" is "-20 meters from the center line"

2) The Translate() and Rotate() functions. Any simple vector math will do here. Please note that if you express your channel profile with reference to 0,0,0 the rotations will be very simple expressions. For more information, please see here: http://en.wikipedia.org/wiki/Rotation_matrix

3) The getAngle() function is a simple atan2(vectorA, vectorB). (for example: http://docs.scipy.org/doc/numpy/reference/generated/numpy.arctan2.html). Please note that in this case, you will be rotating around the Z axis which is the one "sticking out" of your DEM.

4) The orientation of the DEM to the real world. In the above example, we are starting from 0.0 and we will be moving South and then South-East because indexes in a matrix increase top-bottom

Hope this helps a little. Are you dealing with a CFD problem or is this just for visualisation?

CReate 3D elevations over a fixed Grid from a Cross section defined by a string of Points

Setup computational domain

1 Answers