Automatic levee geometries

One of the most boring stuff of levee assessments is the work that has to be done to generate the model for the calculation. It requires the input of the crosssection points, the interpretation of the subsoil and again, the input of all points defining the subsoil.

If something is boring and repetitive then it should probably be automated :-) so this is an article explaining how this can be done using Python.

Ingredients

The things you need are actually quite straightforward. You will need the crosssection as a set of (x,y) points and you will need GEF files (CPT tests for now) as well as the location of the validity of the GEF files. Here is an example;

We have a crosssection and 3 GEF files, the first defining soil section 1, the second section 2 and you get the point. In the script we just expect the right side coordinate of the section and we start the first section at the beginning of the crosssection.

The Python function we use is simple;

def create_geometry(crosssection, gefs, rightlimits)

The crosssection is just a tuple of 2D tuples like

((-10, 10), (0, 10), (3, 1.5).. etc)

In my case I have all GEFs in a local database and I can get them from the database using code like;

for id in (534, 535, 538):
        gefs.append(db.getGEFFromId(id))

The function for the geometry expects GEF objects so you will probably have to write your implementation of that based on my previous article.

The last input parameters is again a simple list or tuple of right side boundaries where the GEF object is valid.

So here is an example of the required input code;

gefs = []
    db = database.Database()
    for id in (534, 535, 538):
        gefs.append(db.getGEFFromId(id))

    create_geometry(crs, gefs, [10, 30, 50])

Now the interesting stuff happens...

The basic Python stuff

Here we go.. we have all input but we need to know the boundary of our geometry. The boundary could be defined by the soillayers found in the CPT or in the crosssection so we check both;

def get_boundary(geometry, crs):
    result = [1e9,-1e9,1e9,-1e9]

    for geo in geometry:
        for layer in geo[2]:
            if layer[1] < result[2]: result[2] = layer[1]
            if layer[1] > result[3]: result[3] = layer[1]

    for p in crs:
        if p[0] < result[0]: result[0] = p[0]
        if p[0] > result[1]: result[1] = p[0]
        if p[1] < result[2]: result[2] = p[1]
        if p[1] > result[3]: result[3] = p[1]

    return tuple(result)

def create_geometry(crosssection, gefs, leftlimits):
    xl = crosssection[0][0]

    geometry = []
    ymin = 1e9
    ymax = -1e9
    for i in range(len(gefs)):
        xr = rightlimits[i]
        layers = []
        ls = gefs[i].toSoil(interval=Y_INTERVAL)
        for l in ls:
            if l[0] > ymax: ymax = l[0]
            if l[1] < ymin: ymin = l[1]
            layers.append(l)
        geometry.append([xl, xr, layers])
        xl = xr

    boundary = get_boundary(geometry, crs)

yes. I agree.. result[1], crosssection[0][0].. not the best code but for a mockup it will do just fine. We will worry about more readable code later ;-)

Note that the..

toSoil()

..function is where I convert the CPT information to soil layers using the simple CUR162 Econe relation. See my previous article for more info. The result is a list of polygons with the name of the soil layer.

The next step is to make sure that our soillayers fit within the boundary. If we find that the first soillayer is lower than the crosssection this might mean that we have no soil between the crosssection and the first layer. We simply change the height of the first soil layer to the maximum of y value and we do the same thing for the last soil layer to minimum of y. Now we are sure that every part of the crosssection is covered.

# correct layers
final_layers = []
for geo in geometry:
    xl, xr = geo[0], geo[1]
    if geo[2][0][0] < boundary[3]: geo[2][0][0] = boundary[3]
    if geo[2][-1][1] > boundary[2]: geo[2][-1][1] = boundary[2]
    for layer in geo[2]:
        final_layers.append((layer[2],((xl, layer[0]),(xr,layer[0]),(xr,layer[1]),(xl,layer[1]))))

Again... [2][-1][1] not the best code but remember.. mockup :-)

We now have rectangles covering our entire geometry.. look at that;

Beautiful but not ready for calculation yet!

The magical Python stuff

The most difficult part is to intersect the rectangles with the crosssection. This crosssection is quite simple but a more complex crosssection and using you own math code would be quite a (fun) challenge. However.. you are not the only one trying to do something with intersections of polygons (which is what we are talking about).. and yes, someone made a c++ library and yes another one made a Python library. Now it's easy;

crosssectionpolygon = update_crs(crosssection, boundary)
intersected_layers = []
for layer in final_layers:
    pc = pyclipper.Pyclipper()
    pc.AddPath(pyclipper.scale_to_clipper(layer[1]), pyclipper.PT_CLIP, True)
    pc.AddPath(pyclipper.scale_to_clipper(crosssectionpolygon), pyclipper.PT_SUBJECT, True)
    sol = pyclipper.scale_from_clipper(pc.Execute(pyclipper.CT_INTERSECTION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD))
    intersected_layers.append((layer[0], tuple([(p[0],p[1]) for p in sol[0]])))

That's all... excellent, let's have a look at the result;

Yes.. that's what we need!

What's next?

The image is in fact made out of a list with entries like this;

('CUR Coarse Sand', ((10.0, 10.0), (0.0, 0.0), (0.0, -2.52), (10.0, -2.52)))

Which is enough information to generate input for calculation software. Unfortunately some commercial software does not support such a polygon approach which is why I am working on my own software that currently is able to make the basic Bishop stability assessment using polygons as soillayers. But that is a whole other story.

I have added the full code here but be aware that it has some dependencies that are not directly available (database as well as the GEF object) so you will have to do something yourself ;-)

Enjoy Python and see you later!

Rob