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?

  • 1
    I'm afraid I don't understand what your input data looks like. Maybe others have the same problem since nobody has turned up with an answer yet.
    – underdark
    Commented Oct 22, 2013 at 17:55
  • I agree could you please be more specific about the input data? You mention various data types, dimensionality, and file formats. Please also consider rephrasing and editing your question for clarity, as you reference the word "string" in different contexts throughout it and the current phrasing is ambiguous. Commented Oct 22, 2013 at 22:28
  • I understand the question, yet can't answer with a solution. As a step forward - here goes my attempt at explaining. Waterways tend to variate in width, depth, etc.. especially when you have lake surfaces in between. Data he has available is basically line data, built from calculated estimates of the water levels (points, connected with a polyline to create a 'surfaceline', or 'string' as the OP described. To visualize in top view which areas are covered in water, he needs to be able to intersect a water surface with a terrain model. Commented May 7, 2014 at 12:58
  • So basically what he needs is to build a 'surface' for the waterline (aka stretch out the line to a surface), convert that into a raster, and compare the DEM with the water surface. Where the DEM is higher, there's no water on the surface, where the DEM is lower then the water level, flooding will occur. Comparison of the two gives an idea of the surface where water will be located, and based on how detailed you can make the comparison, ideally you could create a raster where the value respresents a flood depth. Commented May 7, 2014 at 13:01

1 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])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.