Basic code to add a tree as an excavation in D-Stability using the adjusted geolib version

add_tree_excavation.py

from typing import List, Tuple
import logging
from pydantic import BaseModel
from pathlib import Path
import subprocess
import shutil
from shapely.geometry import Polygon, LineString
from shapely.ops import unary_union, orient

from geolib.models.dstability import DStabilityModel
from geolib.geometry.one import Point

# SETUP DATA LOCATIONS

# Setup the folder to the available calculations
# TODO should point to the sharepoint location
CALCULATIONS_FOLDER = "./data/berekeningen"

# this is a temp directory to save calculation files that are created
# we need to serialize the models before we can calculate them
PATH_TEMP = r"D:\Documents\_TEMP"

# this exe is needed to convert old stix files so geolib can read them
DSTABILITY_BIN_FOLDER = r"D:\Apps\D-GEO Suite\D-Stability 2023.01\bin"

# this allows us to make easy adjustments but the values are static
TREE_WIDTH = 5
TREE_DEPTH = 1

# this script does not use logging but we need to allow some logging
# because we expect problems like missing calculations, faulty
# calculations, trees outside of the calculation limits etc
logging.basicConfig(filename="./ontgrondingskuil.log", level=logging.INFO, filemode="w")


# this is a simple dataclass that is used to store the information
# of the stix file
class StixFile(BaseModel):
    path: str  # path to the file
    dijkvak_id: str  # dijkvak name like A119 based on the filename like A119_0100.stix
    chainage: int  # the chainage (meterering) e.g. 100 based on the filename like A119_0100.stix


def convert_stix(stix_file: str) -> str:
    """Convert the given stix file to the new stix format

    Args:
        stix_file (str): Path to the file

    Returns:
        str: The new path (we add _converted to the original name)
    """
    p = Path(stix_file)
    s = str(Path(p.parent.absolute() / f"{p.stem}_converted.stix"))
    exe = str(Path(DSTABILITY_BIN_FOLDER) / "D-Stability Migration Console.exe")
    subprocess.run([exe, stix_file, s])
    return s


def case_insensitive_glob(filepath: str, fileextension: str) -> List[Path]:
    """Find files with the given extension (e.g. .jpg or .stix) in the given path

    Args:
        filepath (str): The path to search for the files (recursive)
        fileextension (str): The extension to look for (e.g. .jpg, .stix)

    Returns:
        List[Path]: A list of filepaths that are found based on the filter
    """
    p = Path(filepath)
    result = []
    for filename in p.glob("**/*"):
        if str(filename.suffix).lower() == fileextension.lower():
            result.append(filename.absolute())
    return result


def get_stix_files_from_path(path: str) -> List[StixFile]:
    """Find all stix files in the given path and extract extra information
    based on the filename

    Args:
        path (str): Path to the files

    Returns:
        List[StixFile]: A list of StixFile classes with extra information based on the filename
    """
    result = []
    stix_files = case_insensitive_glob(path, ".stix")
    for f in stix_files:
        filename = Path(f).stem

        try:
            args = filename.split("_")
            dtcode = args[0]
            chainage = int(args[1])
            result.append(StixFile(path=str(f), dijkvak_id=dtcode, chainage=chainage))

        except Exception as e:
            # TODO log that this filename cannot be converted to a dijkvak and chainage
            # TODO all files found here need to get a valid filename
            # our team will have to do that
            pass

    return result


def get_surface_from_polygon(polygon: Polygon) -> List[Tuple[float, float]]:
    """Get the surface from a polygon (the highest points from the left to right boundaries)

    Args:
        polygon (Polygon): A shapely polygon

    Raises:
        ValueError: Will be thrown if we cannot find a surface

    Returns:
        List[Tuple[float, float]]: The surface from left to right
    """
    boundary = orient(polygon, sign=-1)
    points = [
        (round(p[0], 3), round(p[1], 3))
        for p in list(zip(*boundary.exterior.coords.xy))[:-1]
    ]

    # get the leftmost point
    left = min([p[0] for p in points])
    right = max([p[0] for p in points])
    topleft_point = sorted([p for p in points if p[0] == left], key=lambda x: x[1])[-1]

    # get the index of this point
    idx = points.index(topleft_point)

    boundary = points[idx:] + points[:idx]
    surface = []
    for p in boundary:
        surface.append(p)
        if p[0] == right:
            return surface

    raise ValueError("Could not determine the surface of the given polygon")


def polygon_from_dstability_model(dm: DStabilityModel) -> Polygon:
    """Convert the dstability model soillayers to one encapsulating polygon,
    Note that we will ONLY convert the first found geometry (scenario 0, stage 0)!

    Args:
        dm (DStabilityModel): The DStabilityModel

    Returns:
        Polygon: A shapely polygon encapsulating all soillayers in the model
    """
    # note we will only use the first geometry in the calculation!
    layers = [
        Polygon([p.X, p.Z] for p in l.Points) for l in dm._get_geometry(0, 0).Layers
    ]
    return unary_union(layers)


def get_z_at(polygon: Polygon, x: float) -> float:
    """Get the highest z coordinate of the given polygon at the given x location

    Args:
        polygon (Polygon): The shapely polygon
        x (float): The location to find the z coordinate

    Raises:
        ValueError: Raised if the x is not within the limits of the polygon

    Returns:
        float: the highest z coordinate of the found intersections
    """
    xx, yy = polygon.exterior.xy
    xmax = max(xx)
    xmin = min(xx)
    ymax = max(yy) + 1.0
    ymin = min(yy) - 1.0

    if x < xmin or x > xmax:
        raise ValueError(
            f"Trying to get a za value outside the x boundary [{xmin},{xmax}]"
        )

    ls = LineString([(x, ymax), (x, ymin)])
    intersections = ls.intersection(polygon)
    return (
        round(max(intersections.xy[1]), 3) + 0.001
    )  # round and add a mm just to be sure


def get_surface_points_between(
    polygon: Polygon, x_start: float, x_end: float
) -> List[Point]:
    """Get the points of the surface of a polygon between the given x coordinates

    Args:
        polygon (Polygon): The shapely polygon
        x_start (float): Start x coordinate
        x_end (float): End x coordinate

    Returns:
        List[Point]: A list of surface points ordered from left to right
    """
    surface = get_surface_from_polygon(polygon)
    return [p for p in surface if p[0] > x_start and p[0] < x_end]


def main():
    # Read in stix files
    stix_files = get_stix_files_from_path(CALCULATIONS_FOLDER)

    # Iterate over all files
    # TODO > this should actually iterate over all trees
    # and per tree find the calculation based on the dijktraject code
    # and chainage
    for stix_file in stix_files:
        # get the original filename
        filename = Path(stix_file.path).stem

        # and create the new one but in the temporary path
        new_filename = str(Path(PATH_TEMP) / f"{filename}.stix")

        # move to temp folder for editing and calculation
        shutil.copyfile(stix_file.path, new_filename)
        stix_file.path = new_filename

        # read the file into the model
        dm = DStabilityModel()
        try:  # old files may create errors, if this happens it probably is an old version
            dm.parse(Path(stix_file.path))
        except ValueError:  # old format, convert and try again
            stix_file.path = convert_stix(stix_file.path)
            try:  # mmm.. nested try except, not really nice but hey, it works!
                dm.parse(Path(stix_file.path))
            except Exception as e:
                pass  # todo, log, ok, that's not the reason.. log the error
        else:  # different error, log
            pass

        # TODO calculate the current SF (safety factor)
        # dm.execute()

        # get the new SF
        # sf_old = dm.get_result(0, 0).FactorOfSafety

        # now the 'difficult part' we don't know the z value at the given x coordinates
        # so we have to do some calculations
        one_big_happy_polygon = polygon_from_dstability_model(dm)

        x_tree = 10  # this should be determined from the database
        x1 = x_tree - TREE_WIDTH / 2.0  # left x coordinate of the 'ontgrondingskuil'
        x2 = x_tree + TREE_WIDTH / 2.0  # right x coordinate of the 'ontgrondingskuil'
        # get the according z coordinates
        z1 = get_z_at(one_big_happy_polygon, x1)
        z2 = get_z_at(one_big_happy_polygon, x2)

        # create a list of points that define the 'ontgrondingskuil'
        # first add the leftmost points
        excavation_points = [Point(x=x1, z=z1), Point(x=x1, z=z1 - TREE_DEPTH)]
        # follow the surface from x1 to x2 and get the intermediate points
        surface_points = get_surface_points_between(one_big_happy_polygon, x1, x2)
        for p in surface_points:
            excavation_points.append(Point(x=p[0], z=p[1] - TREE_DEPTH))
        # add the rightmost points
        excavation_points.append(Point(x=x2, z=z2 - TREE_DEPTH))
        excavation_points.append(Point(x=x2, z=z2))

        # finally.. we can add this excavation using the adjusted geolib version
        dm.add_excavation(label="ontgrondingskuil", points=excavation_points)

        # save this file
        tree_filename = str(Path(PATH_TEMP) / f"{filename}.tree.stix")
        dm.serialize(tree_filename)

        # TODO calculate
        # dm.execute()

        # get the new SF
        # sf_new = dm.get_result(0, 0).FactorOfSafety

        # TODO report SF original and SF new

        break  # remove if you want to run everything although we should actually iterate over the trees :-)


# standard stuff
if __name__ == "__main__":
    main()