add_tree.py

import sys
from pydantic import BaseModel
from pathlib import Path
from typing import Union

# 'Hack' to add the geolib path so Python can find the adjusted geolib code
PATH_TO_ADJUSTED_GEOLIB = r"D:\Documents\Development\Github\GEOLib"
if not PATH_TO_ADJUSTED_GEOLIB in sys.path:
    sys.path.append(PATH_TO_ADJUSTED_GEOLIB)

import geolib as gl
from geolib.models.dstability.loads import TreeLoad, LineLoad
from geolib.geometry.one import Point


class Tree(BaseModel):
    """Definition of a tree for DStability

    Args:
        BaseModel: Inherits from pydantic BaseModel

    Returns:
        Tree: Tree class with the given properties
    """

    height: float  # hoogte
    wind_force: float  # windkracht
    angle_of_distribution: float  # spreidingshoek
    weight: float  # gewicht

    @classmethod
    def from_height(self, height: float) -> "Tree":
        """Generate a tree from the given height. This uses the defaults based on the STOWA publication.
        Note that we are conservative in the assignment of parameters. For example, if a tree if only slightly
        higher than 5m we will use the 10m tree. Also note that the maximum tree is 30m high so anything
        above that height will be cut off at 30m.

        Args:
            height (float): The height of the tree

        Returns:
            Tree: Tree class with the properties according to the STOWA publication
        """
        if height <= 5:
            return Tree(height=5, wind_force=3, angle_of_distribution=30, weight=1)
        elif height <= 10:
            return Tree(height=10, wind_force=6.5, angle_of_distribution=30, weight=5)
        elif height <= 15:
            return Tree(height=15, wind_force=10, angle_of_distribution=30, weight=10)
        elif height <= 20:
            return Tree(height=20, wind_force=13, angle_of_distribution=30, weight=20)
        elif height <= 25:
            return Tree(height=25, wind_force=15, angle_of_distribution=30, weight=40)
        else:
            return Tree(height=30, wind_force=17, angle_of_distribution=30, weight=80)

    @property
    def point_of_force(self) -> float:  # aangrijpingspunt
        return round(self.height * 2.0 / 3.0, 2)


def add_tree(
    stixfile_in: Union[str, Path],
    tree: Tree,
    x: float,
    width_of_rootzone: float = 5,
    stixfile_out: Union[str, Path] = None,
):
    """Add a tree load to a given stix file at the given x coordinate.

    If no stixfile_out parameter is given the stixfile_in will be overwritten else a new stix file will be
    created with the given path and name.

    Args:
        stixfile_in (Union[str, Path]): The path to the file as a string or Path object
        tree (Tree): The tree to add
        x (float): The x coordinate on the geometry in the stix file
        width_of_rootzone (float): The width of the root zone, defaults to 5m
        stixfile_out (Union[str, Path]): The path to the output file, defaults to None which means to overwrite the existing file
    """
    # create the model by reading in the file
    m = gl.DStabilityModel()
    try:
        m.parse(Path(stixfile_in))
    except Exception as e:
        raise ValueError(
            f"Could not open file '{stixfile_in}' for reading. Got error '{e}'"
        )

    # find the z intersections at the given x coordinate
    # note that this uses the adjusted geolib library
    zs = m.datastructure.geometries[0].z_at(x)

    # do we have intersections with the geometry?
    if len(zs) == 0:
        raise ValueError(
            f"There seem to be no intersections with the geometry at the given x coordinate ({x})"
        )

    # get the z coordinate of the top layer and round down so we are sure that it touches the surface
    z = round(zs[0], 2)

    # add a tree load
    # note that this uses the adjusted geolib library
    m.add_load(
        TreeLoad(
            label="tree",
            notes="added by script",
            x=x,
            z=z + tree.point_of_force,
            wind_force=tree.wind_force,
            width_of_root_zone=width_of_rootzone,
            angle_of_distribution=tree.angle_of_distribution,
        )
    )

    # add the lineload which represents the weight of the tree
    m.add_load(
        LineLoad(
            location=Point(x=x, z=z),
            angle=0,
            magnitude=tree.weight,
            angle_of_distribution=tree.angle_of_distribution,
        )
    )

    # overwrite the file?
    if stixfile_out is None:
        stixfile_out = stixfile_in

    # save it
    m.serialize(Path(stixfile_out))


def test():
    """Simple test script"""
    add_tree(
        "./testdata/A147_0230.stix",
        Tree.from_height(10),
        x=0,
        stixfile_out="./testdata/A147_0230_x0.stix",
    )


if __name__ == "__main__":
    test()