Creating water surface from vector line in QGIS

Question

I want to create a raster surface representing the flood level in a creek based on output from a 1D hydraulic model.

I receive a 3D profile of the water surface in a creek line as a 3D Polyline (AutoCad format) from our hydraulic modelers. Each vertex of the polyline represents the coordinates and level of the water surface based on cross-sections used in the modelling.

I want to import the string into QGIS and generate a water surface as a raster for analysis - for example depths, extents, comparison of flood level to building floor levels, etc).

Basically I want to extend a surface horizontally and perpendicular from the string until it intersects a nominated raster surface (existing ground), along the length of the string. I would assume that the elevation varies linearly between each vertex of the polyline (which would be accurate enough for my work).

Is there a way to generate this surface in QGIS?

Answer 1

I think @TimCouwelier’s comments are on the right track.

I used GRASS GIS r.plane to create a raster plane with a given dip (inclination), aspect (azimuth) and one point. You need to provide a list of coordinates with X,Y,Z data called coords and some for loop. This snippet does not provide a full solution but might be a stepping stone into the right direction.

Basically we want to create a raster plane for every set of coordinates along the 3D line. Then as mentioned by @TimCouwelier this result can be compared to the DEM to study where flooding will occur.

#define functions
#Function to calculate the length between the vertices in 3D
def distance_3D(p0, p1):
    return math.sqrt((p1[0] - p0[0])**2 + (p1[1] - p0[1])**2 + (p1[2] - p0[2])**2)

#Function to get the height difference
def delta_Z(p0, p1):
    return (p0[2] - p1[2])

#First calculate the difference in height between the points from a list of vertices
height_dif = delta_Z(coords[0], coords[-1])

#Calculate the length between the vertices in 3D
euclidean_length = distance_3D(coords[0], coords[-1])

#Calculate the dip between two vertices 
dip = math.atan2(height_dif, euclidean_length)

#Calculate azimuth angle
azimuth_angle = abs(180.0/math.pi * math.atan2((coords[0][0] - coords[-1][0]), (coords[0][1] - coords[-1][1])))

#Create plane 
gscript.run_command('r.plane', overwrite=True, output=plane,
    dip=gradient_degree, azimuth=azimuth_angle, easting=coords[-1][0],
    northing=coords[-1][1], elevation=coords[-1][2])