Forestry - TLS

Automate TLS Forestry Analysis - Data Processing

Recap

In the previous tutorial section, we prepared the foundation for automating our TLS-based forestry analysis workflow with Python.

We created a project configuration that centrally defines input data, output directories, filenames, and processing parameters. This allows us to manage and adapt the workflow without modifying the individual processing scripts.

We also developed a data preparation script that creates spatial indices for the input point clouds, establishes a virtual point cloud dataset (.lasvf), and generates the processing units for the buffered and original areas of interest.

With these preparations complete, we are now ready to move on to the actual forestry processing workflow. In this tutorial section, we will use Python to automate the individual processing steps and introduce parallel processing to handle larger TLS datasets more efficiently.

Important

This tutorial section assumes that you have completed the previous section, Automate TLS Forestry Analysis – Data Preparation, and have the project configuration, virtual point cloud dataset, and processing units ready to use.

Create the Processing Scripts

Before creating the processing scripts, let’s briefly revisit the project directory structure established in the previous tutorial section:

project # top-level folder
├── data
├── scripts
|   ├── project_config.py
│   └── 01_data_preparation.py
└── processing

In the scripts folder, create the following four Python scripts:

02_ground_classification.py
03_stem_detection.py
04_tree_segmentation.py
05_tree_models.py

In VS Code, right-click the scripts folder in the Explorer panel and select New File... to create each script.

These scripts will handle the individual processing steps, from ground classification and stem detection to individual tree segmentation and 3D tree reconstruction.

project # top-level folder
├── data
├── scripts
|   ├── project_config.py
|   ├── 01_data_preparation.py
|   ├── 02_ground_classification.py
|   ├── 03_stem_detection.py
|   ├── 04_tree_segmentation.py
│   └── 05_tree_models.py
└── processing

Fill the Ground Classification Script

Open the empty 02_ground_classification.py file in VS Code and copy the following code into it.

Tip

Read the comments in the script to understand the sequence of processing steps performed within the process_unit() function. These comments will help you follow how each processing unit is handled:

# -----------------------------------------------------------------------------
# load processing units

# -----------------------------------------------------------------------------
# load aoi

# -----------------------------------------------------------------------------
# read extent of unit

# -----------------------------------------------------------------------------
# query point cloud for current unit (with overlap)

# -----------------------------------------------------------------------------
# homogenize point cloud

import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import *


@lis_cmd.proc_function_wrapper
def process_unit(unit_index=0):

    # -----------------------------------------------------------------------------
    # load processing units
    processing_units = saga_api.SG_Get_Data_Manager().Add_Shapes(
        cfg.PROC_UNITS_GRDCLASS
    )

    if not processing_units:
        return lis_cmd.Error_Return(
            f"Loading processing_units {cfg.PROC_UNITS_GRDCLASS} failed"
        )

    # -----------------------------------------------------------------------------
    # load aoi
    aoi = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.AOI_BUFFERED)

    if not aoi:
        return lis_cmd.Error_Return(
            f"Loading processing_units {cfg.AOI_BUFFERED} failed"
        )

    # -----------------------------------------------------------------------------
    # read extent of unit
    proc_shape = processing_units.Get_Shape(unit_index)
    xmin = proc_shape.Get_Extent().Get_XMin()
    ymin = proc_shape.Get_Extent().Get_YMin()
    xmax = proc_shape.Get_Extent().Get_XMax()
    ymax = proc_shape.Get_Extent().Get_YMax()

    # -----------------------------------------------------------------------------
    # units are referenced by their lower-left corner coordinate
    unit_name = f"{int(xmin)}_{int(ymin)}"

    x_range = f"{int(xmin)};{int(xmax)}"
    y_range = f"{int(ymin)};{int(ymax)}"

    lis_cmd.Message_Print(f"\nProcessing unit {unit_name} \n")

    # -----------------------------------------------------------------------------
    # select and copy current unit
    processing_units.Select(unit_index)

    proc_unit = saga_api.SG_Get_Data_Manager().Add_Shapes()

    if not shapes_tools.Copy_Selection_to_New_Shapes_Layer(
        INPUT=processing_units, OUTPUT=proc_unit
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    # -----------------------------------------------------------------------------
    # query point cloud for current unit (with overlap)
    pc_out_list = []

    if not lis_virtual.Get_Subset_from_Virtual_LASLAZ(
        AOI_SHP=proc_unit,
        PC_OUT=pc_out_list,
        FILENAME=cfg.LASVF,
        FILEPATH="",
        COPY_ATTR=True,
        CONSTRAIN_QUERY=False,
        ATTR_FIELD="Deviation",
        VALUE_RANGE="0;10",
        FILE_COMPRESSION=True,
        ERROR_IF_EMPTY=False,
        THIN_OUT=False,
        FIELD_TILENAME="<not set>",
        AOI_ADD_OVERLAP=True,
        OVERLAP=cfg.OVERLAP_GRD,
        SKIP_EMPTY_AOIS=False,
        FILENAME_TILE_INFO="",
        METHOD_CLIP="lower and left boundary is inside",
        ONE_PC_PER_POLYGON=False,
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if len(pc_out_list) < 1:
        return lis_cmd.Error_Return(
            f"Processing unit {unit_name} contains no point cloud data, skipping!"
        )

    # add queried point cloud datasets to the data manager
    pointcloud = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()

    if pointcloud.Get_Count() < 1000:
        return lis_cmd.Error_Return(
            f"Processing unit {unit_name} has less than 5 points, skipping ..."
        )

    pointcloud.Get_Projection().Create(cfg.TARGET_COORD)

    # -----------------------------------------------------------------------------
    # homogenize point cloud
    pc_thinned = saga_api.SG_Get_Data_Manager().Add_PointCloud()  # output data object

    if not lis_thinning.Block_Thinning_3D(
        PC_IN=pointcloud,
        PC_OUT=pc_thinned,
        FIELD_THINNING="Reflectance",
        COPY_ATTR=True,
        METHOD_THIN="max(attribute)",
        SPACE_HORI=0.01,
        SPACE_VERT=0.01,
        METHOD_BASELINE="lowest point",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")
    saga_api.SG_Get_Data_Manager().Delete(pointcloud)

    if pc_thinned.Get_Count() < 1000:
        saga_api.SG_Get_Data_Manager().Delete(pc_thinned)
        return lis_cmd.Error_Return(
            f"Tile {unit_name} has less than 1000 points after thinning, skipping ..."
        )

    # -----------------------------------------------------------------------------
    # remove outliers
    pc_filtered = saga_api.SG_Get_Data_Manager().Add_PointCloud()  # output data object

    if not lis_filtering.Remove_Isolated_Points(
        PC_IN=pc_thinned,
        PC_OUT=pc_filtered,
        COPY_ATTR=True,
        SEARCH_RADIUS=cfg.RADIUS_OUTLIER,
        MAX_POINTS=cfg.POINTS_OUTLIER,
        METHOD="remove points",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")
    saga_api.SG_Get_Data_Manager().Delete(pc_thinned)

    if pc_filtered.Get_Count() < 1000:
        saga_api.SG_Get_Data_Manager().Delete(pc_filtered)
        return lis_cmd.Error_Return(
            f"Tile {unit_name} has less than 1000 points after filtering, skipping ..."
        )

    # -----------------------------------------------------------------------------
    # ground classification
    if not lis_classification.Ground_Classification(
        PC_IN=pc_filtered,
        FIELD_CLASS="<not set>",
        FIELD_SEGMENTID="<not set>",
        COPY_ATTR=True,
        SPATIAL_SORT=False,
        ATTR_SUFFIX="",
        INIT_CELLSIZE=2.000000,
        FINAL_CELLSIZE=0.500000,
        AUTOADJUST_CELLSIZE=False,
        FILTER_SEEDS=False,
        TERRAIN_ANGLE=88.000000,
        MIN_EDGE_LENGTH=0.30000,
        MAX_ANGLE=25.000000,
        MAX_DISTANCE=0.800000,
        Z_TOLERANCE=0.100000,
        ATTACH_DZ=True,
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    # -----------------------------------------------------------------------------
    # clip by extent (remove overlap)
    pc_out = saga_api.SG_Get_Data_Manager().Add_PointCloud()

    if not lis_tools.Clip_Point_Cloud_with_Extent(
        PC_IN=pc_filtered,
        PC_OUT=pc_out,
        X_EXTENT=x_range,
        Y_EXTENT=y_range,
        METHOD_CLIP="lower and left boundary is inside",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if pc_out.Get_Count() < 1000:
        return lis_cmd.Error_Return(f"pc_out has not enough points {unit_name}")

    # -----------------------------------------------------------------------------
    # save classified tile (without overlap)
    pc_out.Get_Projection().Create(cfg.TARGET_COORD)
    filename = f"{unit_name}.sg-pts-z"
    pc_out.Save(os.path.join(cfg.CLASS_DIR, filename))

    # -----------------------------------------------------------------------------
    # free memory resources
    saga_api.SG_Get_Data_Manager().Delete()
    return True


# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
if __name__ == "__main__":

    # -----------------------------------------------------------------------------
    # start the logging
    starttime = lis_cmd.Start_Logging(False)

    # -----------------------------------------------------------------------------
    # load the shapefile with the processing units used for processing
    processing_units = saga_api.SG_Get_Data_Manager().Add_Shapes(
        cfg.PROC_UNITS_GRDCLASS
    )

    if not processing_units:
        lis_cmd.Error_Exit(
            f"Loading processing_units {cfg.PROC_UNITS_GRDCLASS} failed",
            starttime,
            True,
        )

    lis_cmd.Message_Print(
        f"Number of units to process: {processing_units.Get_Count()} \n\n"
    )

    # -----------------------------------------------------------------------------
    # run the process_unit function in parallel (each tile in a separate thread)
    # here we will use 8 threads, please adjust to your CPU and available memory
    #
    # the maximum number of threads available with your CPU can be queried with:
    # num_threads = saga_api.SG_OMP_Get_Max_Num_Threads()

    num_threads = 8
    num_processing_units = processing_units.Get_Count()

    lis_cmd.Run_Parallel_Tasks(process_unit, num_processing_units, num_threads)
    # -----------------------------------------------------------------------------

    # -----------------------------------------------------------------------------
    # create a list of all classified point cloud files
    inputfiles = lis_cmd.Create_File_List(
        cfg.CLASS_DIR,
        ".sg-pts;.sg-pts-z",
        os.path.join(cfg.CLASS_DIR, "spcvf_input_file_list.txt"),
        True,
    )

    # -----------------------------------------------------------------------------
    # create virtual point cloud dataset from all classified point cloud files
    if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(
        INPUT_FILE_LIST=inputfiles,
        FILENAME=cfg.SPCVF_CLASS,
        METHOD_PATHS="relative",
        USE_HEADER=True,
        IGNORE_TABLE_STRUCTURE=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # free memory and stop logging
    saga_api.SG_Get_Data_Manager().Delete()
    lis_cmd.Stop_Logging(starttime)

Fill the Stem Detection Script

Open the empty 03_stem_detection.py file in VS Code and copy the following code into it.

import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import *


def process_unit():

    # -----------------------------------------------------------------------------
    # load aoi
    aoi = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.AOI)

    if not aoi:
        return lis_cmd.Error_Return(f"Loading processing_units {cfg.AOI} failed")

    # -----------------------------------------------------------------------------
    # run tree trunk and dbh detection for whole aoi
    if not lis_forestry.Tree_Trunk_and_Diameter_Detection_multitile(
        TILING=aoi,
        SPCVF=cfg.SPCVF_CLASS,
        FIELD_INTENSITY="Reflectance",
        FIELD_DEVIATION="Deviation",
        FIELD_DZ="dz",
        FIELD_TREE_HEIGHT="none",
        COPY_ATTR=True,
        OVERLAP=0.5,
        FILEPATH_SLICE=cfg.SLICES_DIR,
        FILE_COMPRESSION=True,
        FILEPATH_POS=cfg.STEMS_DIR,
        ATTR_SUFFIX="",
        METHOD_STEM_DETECTION="multiple slices",
        SLICE_HEIGHT_RANGE=cfg.SLICE_HEIGHT_RANGE,
        SLICE_STEPPING=cfg.SLICE_STEPPING,
        SLICE_THICKNESS=cfg.SLICE_THICKNESS,
        SLICE_CONNECT_DIST=cfg.SLICE_CONNECT_DIST,
        MIN_INTENSITY=cfg.MIN_INTENSITY,
        MAX_DEVIATION=cfg.MAX_DEVIATION,
        MIN_TREEHEIGHT=cfg.MIN_TREEHEIGHT,
        SEARCH_RADIUS=cfg.SEARCH_RADIUS,
        MIN_POINT_COUNT=cfg.MIN_POINT_COUNT,
        MIN_COMPACTNESS=cfg.MIN_COMPACTNESS,
        STEM_DIAMETER_RANGE=cfg.STEM_DIAMETER_RANGE,
        CIRCLE_FIT_TOL=cfg.CIRCLE_FIT_TOL,
        POINTS_ON_CIRCLE=cfg.POINTS_ON_CIRCLE,
        GOODNESS_OF_FIT=cfg.GOODNESS_OF_FIT,
        CIRCLE_COMPLETENESS=cfg.CIRCLE_COMPLETENESS,
        MIN_DETECTIONS=cfg.MIN_DETECTIONS,
        PERCENTILE_RADIUS=20.000000,
        APPLY_CORR_FUNC=False,
    ):
        lis_cmd.Error_Return("Failed to execute dbh detection tool on tile")

    # -----------------------------------------------------------------------------
    # free memory resources
    saga_api.SG_Get_Data_Manager().Delete()
    return True


# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
if __name__ == "__main__":

    # -----------------------------------------------------------------------------
    # start the logging
    starttime = lis_cmd.Start_Logging(False)

    # -----------------------------------------------------------------------------
    # call the processing function
    result = process_unit()

    if not result:
        lis_cmd.Error_Exit("Processing task failed", starttime)

    # -----------------------------------------------------------------------------
    # create a list of all point cloud slice files
    inputfiles = lis_cmd.Create_File_List(
        cfg.SLICES_DIR,
        ".sg-pts;.sg-pts-z",
        os.path.join(cfg.SLICES_DIR, "spcvf_input_file_list.txt"),
        True,
    )

    # -----------------------------------------------------------------------------
    # create virtual point cloud from all point cloud slice files
    if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(
        INPUT_FILE_LIST=inputfiles,
        FILENAME=cfg.SPCVF_SLICES,
        METHOD_PATHS="relative",
        USE_HEADER=True,
        IGNORE_TABLE_STRUCTURE=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # load stem positions
    stem_pos = saga_api.SG_Get_Data_Manager().Add_Shapes(
        os.path.join(cfg.STEMS_DIR, "tree_positions.shp")
    )

    if not stem_pos:
        lis_cmd.Error_Exit("Loading stem positions failed", starttime, True)

    # -----------------------------------------------------------------------------
    # load reference stem positions
    # optional step in case you already have georeferenced tree positions from a catastre
    # reference positions could also have reference tree metrics such as dbh and tree height
    reference_pos = saga_api.SG_Get_Data_Manager().Add_Shapes(
        os.path.join(cfg.REFERENCE_TREES)
    )

    if not reference_pos:
        lis_cmd.Error_Exit("Loading reference stem positions failed", starttime, True)

    # -----------------------------------------------------------------------------
    # join the provided reference tree ids to the derived detections (raw tree ids are processing dependent (processing ids))
    if not lis_tools.Join_Attributes_by_NN_Search_Shapes(
        SHAPES_IN=stem_pos,
        SHAPES_JOIN=reference_pos,
        RESULT=result,
        FIELD_MAX_DISTANCE="<not set>",
        FIELDS="ID",
        ATTR_SUFFIX="Reference",
        DIMENSION="2D search",
        MAX_DISTANCE=0.200000,
        METHOD_UNDEFINED="NoData",
        ADD_JOIN_FLAG=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # update stem pos with reference tree id
    stem_pos.Save(os.path.join(cfg.STEMS_DIR, "tree_positions.shp"))

    # -----------------------------------------------------------------------------
    # load aoi
    aoi = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.AOI)

    if not aoi:
        lis_cmd.Error_Exit("Loading aoi failed", starttime, True)

    # -----------------------------------------------------------------------------
    # generate dbh slice grid for the whole area of interest
    grids_list = []

    if not lis_virtual.Get_Grid_from_Virtual_Point_Cloud(
        AOI_SHP=aoi,
        GRID_OUT=grids_list,
        FILENAME=cfg.SPCVF_SLICES,
        FILEPATH="",
        ATTR_FIELD_GRID="Reflectance",
        CELL_SIZE=0.010000,
        GRID_SYSTEM_FIT="cells",
        METHOD="maximum",
        CONSTRAIN_QUERY=False,
        FIELD_TILENAME="<not set>",
        AOI_XRANGE="0; 0",
        AOI_YRANGE="0; 0",
        AOI_ADD_OVERLAP=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    if len(grids_list) < 1:
        lis_cmd.Error_Exit("slices are empty", starttime, True)

    # add the output grid dataset to the data manager
    grid = saga_api.SG_Get_Data_Manager().Add(grids_list[0]).asGrid()

    # -----------------------------------------------------------------------------
    # save the slice grid
    grid.Get_Projection().Create(cfg.TARGET_COORD)
    grid.Set_Name(f"{cfg.PLOT_NAME}_slices_grid")
    grid.Save(os.path.join(cfg.GRIDS_DIR, f"{cfg.PLOT_NAME}_slices_grid.tif"))

    # -----------------------------------------------------------------------------
    # free memory and stop logging
    saga_api.SG_Get_Data_Manager().Delete()
    lis_cmd.Stop_Logging(starttime)

Fill the Tree Segmentation Script

Open the empty 04_tree_segmentation.py file in VS Code and copy the following code into it.

import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import *


@lis_cmd.proc_function_wrapper
def process_unit(unit_index=0):

    # -----------------------------------------------------------------------------
    # load processing units
    processing_units = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.PROC_UNITS)

    if not processing_units:
        return lis_cmd.Error_Return(f"Loading processing_units {cfg.PROC_UNITS} failed")

    globalMinX = processing_units.Get_Extent().Get_XMin()
    globalMinY = processing_units.Get_Extent().Get_YMin()

    # -----------------------------------------------------------------------------
    # load aoi (without buffer)
    aoi = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.AOI)

    if not aoi:
        return lis_cmd.Error_Return(f"Loading aoi {cfg.AOI} failed")

    # -----------------------------------------------------------------------------
    # load stem positions
    stem_pos = saga_api.SG_Get_Data_Manager().Add_Shapes(
        os.path.join(cfg.STEMS_DIR, "tree_positions.shp")
    )

    if not stem_pos:
        return lis_cmd.Error_Return(f"Loading tree positions failed")

    stem_pos_selected = (
        saga_api.SG_Get_Data_Manager().Add_Shapes()
    )  # output data object

    if not shapes_tools.Select_by_Attributes_Numerical_Expression(
        SHAPES=stem_pos,
        COPY=stem_pos_selected,
        FIELD="<not set>",
        EXPRESSION="[ID_Reference] > 0",
        INVERSE=False,
        USE_NODATA=False,
        METHOD="new selection",
        POSTJOB="copy",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    # -----------------------------------------------------------------------------
    # select and copy current processing unit
    processing_units.Select(unit_index)

    proc_unit = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object

    if not shapes_tools.Copy_Selection_to_New_Shapes_Layer(
        INPUT=processing_units, OUTPUT=proc_unit
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    xmin = proc_unit.Get_Extent().Get_XMin()
    ymin = proc_unit.Get_Extent().Get_YMin()

    unit_name = f"{int(xmin)}_{int(ymin)}"

    # -----------------------------------------------------------------------------
    # add buffer to the processing unit
    buffer = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object

    if not shapes_tools.Shapes_Buffer(
        SHAPES=proc_unit,
        BUFFER=buffer,
        DIST_FIELD="<not set>",
        DIST_FIELD_DEFAULT=cfg.OVERLAP,
        DISSOLVE=True,
        NZONES=1,
        DARC=5.000000,
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    # -----------------------------------------------------------------------------
    # query point cloud for buffered processing unit
    pc_out_list = []

    if not lis_virtual.Get_Subset_from_Virtual_Point_Cloud(
        AOI_SHP=buffer,
        PC_OUT=pc_out_list,
        FILENAME=cfg.SPCVF_CLASS,
        FILEPATH="",
        FILE_COMPRESSION=True,
        COPY_ATTR=True,
        CONSTRAIN_QUERY=False,
        THIN_OUT=False,
        ERROR_IF_EMPTY=False,
        FIELD_TILENAME="<not set>",
        AOI_ADD_OVERLAP=False,
        OVERLAP=0,
        METHOD_CLIP="lower and left boundary is inside",
        ONE_PC_PER_POLYGON=False,
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if len(pc_out_list) < 1:
        return lis_cmd.Error_Return("Failed to execute tool - output list empty")

    pc_buff = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()

    if pc_buff.Get_Count() < 5:
        return lis_cmd.Error_Return("Failed to execute tool - not enough points")

    # -----------------------------------------------------------------------------
    # single tree segmentation
    tree_positions = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object
    boundaries = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object

    if not (
        lis_forestry.Single_Tree_Derivation_3D(
            PC_IN=pc_buff,
            TREE_POS_IN=stem_pos_selected,
            TREE_POS_OUT=tree_positions,
            BOUNDARIES=boundaries,
            FIELD_INTENSITY="Reflectance",
            FIELD_DEVIATION="Deviation",
            FIELD_CLASS="grdclass",
            FIELD_DZ="dz",
            COPY_ATTR=True,
            EXCLUDE_CLASSES="2;6;9",
            FIELD_TREEID="ID_Reference",  # reference treeid
            FIELD_TRUNK_DIAMETER="diameter_mean",
            FIELD_SLICE_HEIGHT="slice_height",
            FIELD_BLOCK_DIST="<not set>",
            ATTR_SUFFIX="",
            SEED_STEPPING=cfg.SEED_STEPPING,
            WEIGHT=cfg.WEIGHT,
            MAX_GAP_SIZE=cfg.MAX_GAP_SIZE,
            FOLIAGE_DISTANCE=cfg.FOLIAGE_DISTANCE,
            MIN_INTENSITY=cfg.MIN_INTENSITY_SEG,
            MAX_DEVIATION=cfg.MAX_DEVIATION_SEG,
            MIN_PT_DENSITY=cfg.MIN_PT_DENSITY,
            MIN_TREE_HEIGHT=cfg.MIN_TREE_HEIGHT,
            MAX_LOWVEG_HEIGHT=cfg.MAX_LOWVEG_HEIGHT,
            CROWN_METRICS_RES=cfg.CROWN_METRICS_RES,
            NUM_ITERATIONS=cfg.NUM_ITERATIONS,
            ALL_TREES=cfg.ALL_TREES,
            START_COST_SYNTH=0.000000,
            SPECIFY_ORIGIN=False,
            TREEID_OFFSET=0,
            LOWERLEFT_X=globalMinX,
            LOWERLEFT_Y=globalMinY,
        )
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if tree_positions.Get_Count() < 1:
        return lis_cmd.Error_Return(
            "Failed to execute tool - not enough tree positions"
        )

    # -----------------------------------------------------------------------------
    # drop unneccesary attributes
    if not lis_tools.Drop_Point_Cloud_Attributes(
        PC_IN=pc_buff, FIELDS="distdiff,dist_to_seed,voxel_pt_cnt,pathdist"
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    # -----------------------------------------------------------------------------
    # clip by processing unit border (without buffer)
    pc_tile = saga_api.SG_Get_Data_Manager().Add_PointCloud()  # output data object
    tree_positions_tile = (
        saga_api.SG_Get_Data_Manager().Add_Shapes()
    )  # output data object

    if not lis_forestry.Tree_Segment_Aware_Overlap_Removal(
        PC_IN=pc_buff,
        STEM_POSITIONS=tree_positions,
        SELECTION_POLYGON=proc_unit,
        PC_OUT=pc_tile,
        STEM_POSITIONS_OUT=tree_positions_tile,
        FIELD_TREEID="treeid",
        FIELD_TREEID_STEM="treeid",
        FIELD_X_STEM="<not set>",
        FIELD_Y_STEM="<not set>",
        FIELD_TREEID_SKEL="<no attributes>",
        FIELD_TREEID_PIPE="<no attributes>",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if tree_positions_tile.Get_Count() < 1:
        return lis_cmd.Error_Return("Not enough trees after clipping")

    if pc_tile.Get_Count() < 5:
        return lis_cmd.Error_Return("Not enough points after clipping")

    # -----------------------------------------------------------------------------
    # clip by aoi
    pc_out = saga_api.SG_Get_Data_Manager().Add_PointCloud()  # output data object
    tree_positions_out = (
        saga_api.SG_Get_Data_Manager().Add_Shapes()
    )  # output data object

    if not lis_forestry.Tree_Segment_Aware_Overlap_Removal(
        PC_IN=pc_tile,
        STEM_POSITIONS=tree_positions_tile,
        SELECTION_POLYGON=aoi,
        PC_OUT=pc_out,
        STEM_POSITIONS_OUT=tree_positions_out,
        FIELD_TREEID="treeid",
        FIELD_TREEID_STEM="treeid",
        FIELD_X_STEM="<not set>",
        FIELD_Y_STEM="<not set>",
        FIELD_TREEID_SKEL="<no attributes>",
        FIELD_TREEID_PIPE="<no attributes>",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if tree_positions_out.Get_Count() < 1:
        return lis_cmd.Error_Return("Not enough trees after clipping")

    if pc_out.Get_Count() < 5:
        return lis_cmd.Error_Return("Not enough points after clipping")

    # -----------------------------------------------------------------------------
    # select crown boundaries corresponding to current processing unit
    boundaries_out = saga_api.SG_Get_Data_Manager().Add_Shapes()

    if not shapes_tools.Select_by_Location(
        SHAPES=boundaries,
        COPY=boundaries_out,
        LOCATIONS=tree_positions_out,
        CONDITION="intersect",
        METHOD="new selection",
        POSTJOB="copy",
    ):
        return lis_cmd.Error_Return("Failed to execute tool")

    if boundaries_out.Get_Count() < 1:  ## there is nothing to keep
        return lis_cmd.Error_Return("Not enough tree boundaries after clipping")

    # -----------------------------------------------------------------------------
    # save boundaries
    boundaries_name = f"{unit_name}_boundaries.shp"
    boundaries_out.Get_Projection().Create(cfg.TARGET_COORD)
    boundaries_out.Save(os.path.join(cfg.OUTLINES_DIR, boundaries_name))

    # -----------------------------------------------------------------------------
    # save tree positions
    positions_name = f"{unit_name}_positions.shp"
    tree_positions_out.Get_Projection().Create(cfg.TARGET_COORD)
    tree_positions_out.Save(os.path.join(cfg.TREE_POS_DIR, positions_name))

    # -----------------------------------------------------------------------------
    # save segmented point cloud tiles
    pc_name = f"{unit_name}_segments.sg-pts-z"
    pc_out.Get_Projection().Create(cfg.TARGET_COORD)
    pc_out.Save(os.path.join(cfg.SEGMENTS_DIR, pc_name))

    # -----------------------------------------------------------------------------
    # free memory ressources
    saga_api.SG_Get_Data_Manager().Delete()
    return True


# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
if __name__ == "__main__":

    # -----------------------------------------------------------------------------
    # start the logging
    starttime = lis_cmd.Start_Logging(False)

    # -----------------------------------------------------------------------------
    # load processing units
    processing_units = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.PROC_UNITS)

    if not processing_units:
        lis_cmd.Error_Exit(
            f"Loading processing_units {cfg.PROC_UNITS} failed",
            starttime,
            True,
        )

    lis_cmd.Message_Print(
        f"Number of units to process: {processing_units.Get_Count()} \n\n"
    )

    # -----------------------------------------------------------------------------
    # run the process_unit function in parallel (each tile in a separate thread)
    # here we will use 4 threads, please adjust to your CPU and available memory
    #
    # the maximum number of threads available with your CPU can be queried with:
    # num_threads = saga_api.SG_OMP_Get_Max_Num_Threads()

    num_threads = 4
    num_processing_units = processing_units.Get_Count()

    lis_cmd.Run_Parallel_Tasks(process_unit, num_processing_units, num_threads)
    # -----------------------------------------------------------------------------

    # -----------------------------------------------------------------------------
    # create a list of all segmented point cloud files
    inputfiles = lis_cmd.Create_File_List(
        cfg.SEGMENTS_DIR,
        ".sg-pts;.sg-pts-z",
        os.path.join(cfg.SEGMENTS_DIR, "spcvf_input_file_list.txt"),
        True,
    )

    # -----------------------------------------------------------------------------
    # create virtual point cloud from segmented point cloud files
    if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(
        INPUT_FILE_LIST=inputfiles,
        FILENAME=cfg.SPCVF_SEGMENTS,
        METHOD_PATHS="relative",
        USE_HEADER=True,
        IGNORE_TABLE_STRUCTURE=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # create a list of all tree position files
    pathlist = lis_cmd.Create_File_List(
        cfg.TREE_POS_DIR,
        ".shp;.SHP",
        os.path.join(cfg.TREE_POS_DIR, "shp_input_file_list.txt"),
        False,
    )

    # -----------------------------------------------------------------------------
    # load all tree position files into a python list object
    merge_list = []
    for file in lis_cmd.Get_File_List_As_List(pathlist):

        positions = saga_api.SG_Get_Data_Manager().Add_Shapes(file)

        if not positions or positions.Get_Count() < 1:
            lis_cmd.Message_Print(f"Loading {file} failed")
            continue

        merge_list.append(positions)

    # -----------------------------------------------------------------------------
    # merge all loaded tree position datasets into one
    merged_positions = saga_api.SG_Get_Data_Manager().Add_Shapes()

    if not shapes_tools.Merge_Layers(
        INPUT=merge_list,
        MERGED=merged_positions,
        SRCINFO=False,
        MATCH=True,
        DELETE=True,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # save merged tree positions
    merged_positions.Save(os.path.join(cfg.SEGMENTATION_DIR, "tree_positions.shp"))

    # -----------------------------------------------------------------------------
    # create a list of all available crown outline files
    pathlist = lis_cmd.Create_File_List(
        cfg.OUTLINES_DIR,
        ".shp;.SHP",
        os.path.join(cfg.OUTLINES_DIR, "shp_input_file_list.txt"),
        False,
    )

    # -----------------------------------------------------------------------------
    # load all crown outline files into a python list object
    merge_list = []  # the list holding all input objects to be merged
    for file in lis_cmd.Get_File_List_As_List(pathlist):

        outlines = saga_api.SG_Get_Data_Manager().Add_Shapes(file)

        if not outlines or outlines.Get_Count() < 1:
            lis_cmd.Message_Print(f"Loading {file} failed")
            continue
        merge_list.append(outlines)

    # -----------------------------------------------------------------------------
    # merge all loaded outline datasets into one
    merged_outlines = saga_api.SG_Get_Data_Manager().Add_Shapes()
    if not shapes_tools.Merge_Layers(
        INPUT=merge_list, MERGED=merged_outlines, SRCINFO=False, MATCH=True, DELETE=True
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # save merged crown outlines
    merged_outlines.Save(os.path.join(cfg.SEGMENTATION_DIR, "tree_outlines.shp"))

    # -----------------------------------------------------------------------------
    # load aoi
    aoi = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.AOI_BUFFERED)

    if not aoi:
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # generate chm grid for the whole area of interest
    grids_list = []

    if not lis_virtual.Get_Grid_from_Virtual_Point_Cloud(
        AOI_SHP=aoi,
        GRID_OUT=grids_list,
        FILENAME=cfg.SPCVF_SEGMENTS,
        FILEPATH="",
        ATTR_FIELD_GRID="dz",
        CELL_SIZE=0.1,
        GRID_SYSTEM_FIT="cells",
        METHOD="maximum",
        CONSTRAIN_QUERY=False,
        FIELD_TILENAME="<not set>",
        AOI_XRANGE="0; 0",
        AOI_YRANGE="0; 0",
        AOI_ADD_OVERLAP=False,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    if len(grids_list) < 1:
        lis_cmd.Error_Exit("slices are empty", starttime, True)

    # add the output grid dataset to the data manager
    chm = saga_api.SG_Get_Data_Manager().Add(grids_list[0]).asGrid()

    # -----------------------------------------------------------------------------
    # save chm
    chm.Get_Projection().Create(cfg.TARGET_COORD)
    chm.Set_Name(f"{cfg.PLOT_NAME}_chm")
    chm.Save(os.path.join(cfg.GRIDS_DIR, f"{cfg.PLOT_NAME}_chm.tif"))

    # -----------------------------------------------------------------------------
    # free memory and stop logging
    saga_api.SG_Get_Data_Manager().Delete()
    lis_cmd.Stop_Logging(starttime)

Fill the Tree Models Script

Open the empty 05_tree_models.py file in VS Code and copy the following code into it.

import os
import csv
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import *


def recode_fields(inpath):
    outpath = os.path.join(os.path.dirname(inpath), "tmp.csv")

    with open(inpath, "r", newline="", encoding="utf-8") as infile, open(
        outpath, "w", newline="", encoding="utf-8"
    ) as outfile:

        reader = csv.DictReader(infile)
        writer = csv.DictWriter(outfile, fieldnames=reader.fieldnames)

        writer.writeheader()

        for row in reader:

            # ----------------------------------------
            # replace nodata
            for field, value in row.items():
                try:
                    if float(value) == -1.0:
                        row[field] = ""
                except ValueError:
                    pass

            writer.writerow(row)


@lis_cmd.proc_function_wrapper
def process_unit(unit_index=0):

    # -----------------------------------------------------------------------------
    # Create a list of all segmented tiles.
    # Tiles already include complex borders to avoid cutting tree segments.
    input_filelist = lis_cmd.Get_File_List_As_List(
        os.path.join(cfg.SEGMENTS_DIR, "spcvf_input_file_list.txt")
    )

    # -----------------------------------------------------------------------------
    # get current tile by unit_index.
    file = input_filelist[unit_index]

    basename = os.path.basename(file)
    unit_name, unitext = os.path.splitext(basename)

    unit_name = unit_name.replace("_segments", "")

    # --------------------------------------------------------------
    # load point cloud of processing unit
    pc = saga_api.SG_Get_Data_Manager().Add_PointCloud(file)

    if not pc:
        return lis_cmd.Error_Return(f"Loading point cloud {file} failed")

    if pc.Get_Count() < 1000:
        return lis_cmd.Error_Return("point cloud has not enough points, skipping... ")

    # --------------------------------------------------------------
    # load tree positions of processing unit
    positions_name = f"{unit_name}_positions.shp"

    tree_pos = saga_api.SG_Get_Data_Manager().Add_Shapes(
        os.path.join(os.path.join(cfg.TREE_POS_DIR, positions_name))
    )

    if not tree_pos:
        return lis_cmd.Error_Return(f"Loading processing_units tree positions failed")

    # --------------------------------------------------------------
    # compute tree shape metrics
    tree_metrics = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object
    shp_skeletons = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object
    shp_pipes = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object
    shp_branches = saga_api.SG_Get_Data_Manager().Add_Shapes()  # output data object

    if not lis_forestry.Tree_Shape_Metrics_3D(
        PC_IN=pc,
        TREE_POS_IN=tree_pos,
        TREE_METRICS=tree_metrics,
        SHP_SKELETON=shp_skeletons,
        SHP_SKELETON_SIMPLE=shp_branches,
        SHP_PIPE=shp_pipes,
        FIELD_INTENSITY="Reflectance",
        FIELD_DEVIATION="Deviation",
        FIELD_TREEID="treeid",
        COPY_ATTR=True,
        FIELD_TREEID_TREE_POS="treeid",
        FIELD_TRUNK_DIAMETER="trunk_diameter",
        FIELD_SLICE_HEIGHT="slice_height",
        DIAMETER_FACTOR=cfg.DIAMETER_FACTOR,
        ATTR_SUFFIX="",
        MIN_INTENSITY=cfg.MIN_INTENSITY_METRICS,
        MAX_DEVIATION=cfg.MAX_DEVIATION_METRICS,
        INT_VALUE_RANGE=cfg.INT_VALUE_RANGE,
        INT_SCALE_FACTOR=cfg.INT_SCALE_FACTOR,
        SEARCH_RADIUS=cfg.SEARCH_RADIUS_METRICS,
        BRIDGE_DIST=cfg.BRIDGE_DIST,
        ITERATIONS=cfg.ITERATIONS,
        NODE_STEPPING=cfg.NODE_STEPPING,
        BRANCHING_DIST=cfg.BRANCHING_DIST,
        NUMBER_PRISMFACES=cfg.NUMBER_PRISMFACES,
        MIN_CHILDREN=cfg.MIN_CHILDREN,
        MIN_DIAMETER=cfg.MIN_DIAMETER,
        FORK_RATIO=cfg.FORK_RATIO,
        MIN_FORK_DIAMETER=cfg.MIN_FORK_DIAMETER,
        SIGMA_SMOOTHING=cfg.SIGMA_SMOOTHING,
        DIRECTION_DETECTION_LENGTH=cfg.DIRECTION_DETECTION_LENGTH,
        DBH_REFERENCE_HEIGHT=cfg.DBH_REFERENCE_HEIGHT,
        ADJUSTMENT_MODE=cfg.ADJUSTMENT_MODE,
        ADJUSTMENT_STEPPING=cfg.ADJUSTMENT_STEPPING,
        ADJUSTMENT_START=cfg.ADJUSTMENT_START,
        TRUNK_BASE_HEIGHT=cfg.TRUNK_BASE_HEIGHT,
        CROSS_SECTION_MODE=cfg.CROSS_SECTION_MODE,
        CIRCULAR_CALC_METHOD=cfg.CIRCULAR_CALC_METHOD,
        REFINE_BRANCH_NODES=cfg.REFINE_BRANCH_NODES,
        SEARCH_RADIUS_COMPACTNESS=cfg.SEARCH_RADIUS_COMPACTNESS,
        MIN_COMPACTNESS=cfg.MIN_COMPACTNESS_METRICS,
        CIRCLE_FIT_TOL=cfg.CIRCLE_FIT_TOL_METRICS,
        CIRCLE_FIT_RESOLUTION=cfg.CIRCLE_FIT_RESOLUTION,
        DIAMETER_CORR_FACTOR=cfg.DIAMETER_CORR_FACTOR,
        GOODNESS_OF_FIT=cfg.GOODNESS_OF_FIT_METRICS,
        CIRCLE_COMPLETENESS=cfg.CIRCLE_COMPLETENESS_METRICS,
        TAPER_TOLERANCE=cfg.TAPER_TOLERANCE,
        THINNING=True,
        LONG_LOG_LENGTH=6.000000,
        HQ_STRAIGHT_TOL=8.000000,
        MQ_STRAIGHT_TOL_LONG=5.000000,
        SHORT_LOG_LENGTH=2.000000,
        MQ_STRAIGHT_TOL_SHORT=5.000000,
        AQ_STRAIGHT_TOL=3.000000,
    ):
        return lis_cmd.Error_Return("Failed to execute tool on tile")

    # --------------------------------------------------------------
    # save tree metrics
    if tree_metrics.Get_Count() > 0:
        metrics_name = f"{unit_name}_metrics.shp"

        tree_metrics.Get_Projection().Create(cfg.TARGET_COORD)
        tree_metrics.Save(os.path.join(cfg.TREE_METRICS_DIR, metrics_name))

    # --------------------------------------------------------------
    # save tree skeletons
    if shp_skeletons.Get_Count() > 0:
        skeletons_name = f"{unit_name}_skeletons.shp"

        shp_skeletons.Get_Projection().Create(cfg.TARGET_COORD)
        shp_skeletons.Save(os.path.join(cfg.SKELETONS_DIR, skeletons_name))

    # --------------------------------------------------------------
    # save branches
    if shp_branches.Get_Count() > 0:
        branches_name = f"{unit_name}_branches.shp"

        shp_branches.Get_Projection().Create(cfg.TARGET_COORD)
        shp_branches.Save(os.path.join(cfg.BRANCH_AXES_DIR, branches_name))

    # --------------------------------------------------------------
    # save pipe models
    if shp_pipes.Get_Count() > 0:

        # -----------------------------------------------------------------------------
        # extract higher geometries of pipe model
        query = f"[Hierarchy] < {cfg.MAX_HIERARCHY_REPORT+1}"
        shp_pipes_main = (
            saga_api.SG_Get_Data_Manager().Add_Shapes()
        )  # output data object

        if not shapes_tools.Select_by_Attributes_Numerical_Expression(
            SHAPES=shp_pipes,
            COPY=shp_pipes_main,
            FIELD="<not set>",
            EXPRESSION=query,
            INVERSE=False,
            USE_NODATA=False,
            METHOD="new selection",
            POSTJOB="copy",
        ):
            return lis_cmd.Error_Return("Failed to execute tool")

        if shp_pipes_main.Get_Count() < 1:
            return lis_cmd.Error_Return(f"single tree pipemodel has no entries")
        pipes_name = f"{unit_name}_pipes.shp"
        shp_pipes_main.Get_Projection().Create(cfg.TARGET_COORD)
        shp_pipes_main.Save(os.path.join(cfg.PIPES_DIR, pipes_name))

    # -----------------------------------------------------------------------------
    # free memory ressources
    saga_api.SG_Get_Data_Manager().Delete()
    return True


# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
if __name__ == "__main__":

    # -----------------------------------------------------------------------------
    # start the logging
    starttime = lis_cmd.Start_Logging(False)

    # -----------------------------------------------------------------------------
    # create a Python list of all segmented tiles
    file_list = lis_cmd.Get_File_List_As_List(
        os.path.join(cfg.SEGMENTS_DIR, "spcvf_input_file_list.txt")
    )

    # -----------------------------------------------------------------------------
    # run the process_unit function in parallel (each tile in a separate thread)
    # here we will use 4 threads, please adjust to your CPU and available memory
    #
    # the maximum number of threads available with your CPU can be queried with:
    # num_threads = saga_api.SG_OMP_Get_Max_Num_Threads()

    num_threads = 4
    num_processing_units = len(file_list)

    lis_cmd.Run_Parallel_Tasks(process_unit, num_processing_units, num_threads)
    # -----------------------------------------------------------------------------

    # -----------------------------------------------------------------------------
    # create a list of all tree metrics files
    pathlist = lis_cmd.Create_File_List(
        cfg.TREE_METRICS_DIR,
        ".shp;.SHP",
        os.path.join(cfg.TREE_METRICS_DIR, "shp_input_file_list.txt"),
        False,
    )

    # -----------------------------------------------------------------------------
    # load all tree metric files into a Python list object
    merge_list = []
    for file in lis_cmd.Get_File_List_As_List(pathlist):

        metrics = saga_api.SG_Get_Data_Manager().Add_Shapes(file)

        if not metrics or metrics.Get_Count() < 1:
            lis_cmd.Message_Print(f"Loading {file} failed")
            continue
        merge_list.append(metrics)

    # -----------------------------------------------------------------------------
    # merge all loaded tree metric datasets into one
    merged_metrics = saga_api.SG_Get_Data_Manager().Add_Shapes()
    if not shapes_tools.Merge_Layers(
        INPUT=merge_list, MERGED=merged_metrics, SRCINFO=False, MATCH=True, DELETE=True
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # save merged tree metrics
    merged_metrics.Save(os.path.join(cfg.METRICS_DIR, "tree_metrics.shp"))

    # -----------------------------------------------------------------------------
    # create a list all branch axes files
    pathlist = lis_cmd.Create_File_List(
        cfg.BRANCH_AXES_DIR,
        ".shp;.SHP",
        os.path.join(cfg.BRANCH_AXES_DIR, "shp_input_file_list.txt"),
        False,
    )

    # -----------------------------------------------------------------------------
    # load all branch axes files into a Python list object
    merge_list = []
    for file in lis_cmd.Get_File_List_As_List(pathlist):

        branches = saga_api.SG_Get_Data_Manager().Add_Shapes(file)

        if not branches.is_Valid() or branches.Get_Count() < 1:
            lis_cmd.Message_Print(f"Loading {file} failed")
            continue
        merge_list.append(branches)

    # -----------------------------------------------------------------------------
    # merge all loaded branch axes datasets into one
    merged_branches = saga_api.SG_Get_Data_Manager().Add_Shapes()
    if not shapes_tools.Merge_Layers(
        INPUT=merge_list, MERGED=merged_branches, SRCINFO=False, MATCH=True, DELETE=True
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # save merged branches
    merged_branches.Save(os.path.join(cfg.METRICS_DIR, "tree_branches.shp"))

    # -----------------------------------------------------------------------------
    # distill tree information csv tables
    tree_tab = saga_api.SG_Get_Data_Manager().Add_Table()  # output data object
    stem_tab = saga_api.SG_Get_Data_Manager().Add_Table()  # output data object

    if not lis_forestry.Extract_Tree_Information_3D(
        TREE_METRICS=merged_metrics,
        BRANCHES=merged_branches,
        TREE_TAB=tree_tab,
        STEM_TAB=stem_tab,
        FIELD_TREEID="treeid",
        FIELD_DBH="DBH",
        FIELD_TREEHEIGHT="tree_height",
        FIELD_CROWNDIAMETER="crown_diameter",
        FIELD_CROWNAREA="crown_area",
        FIELD_CROWNBOTTOM="CrownBottom",
        FIELD_CROWN_MAXRADIUS="crown_max_rad",
        POPULATION_NAME=cfg.POPULATION_NAME,
        STRATUM_NAME=cfg.STRATUM_NAME,
        PLOT_NAME=cfg.PLOT_NAME,
        DBH_HEIGHT=1.3,
        ANALYZE_BRANCHES=False,
        CALC_BIOMASS=False,
        ADD_CROWN_INFO=True,
    ):
        lis_cmd.Error_Exit("", starttime, True)

    # -----------------------------------------------------------------------------
    # save tables
    tree_tab.Save(os.path.join(cfg.TABLES_DIR, "tree_info.csv"))
    stem_tab.Save(os.path.join(cfg.TABLES_DIR, "stem_description.csv"))
    recode_fields(os.path.join(cfg.TABLES_DIR, "tree_info.csv"))
    recode_fields(os.path.join(cfg.TABLES_DIR, "stem_description.csv"))

    # -----------------------------------------------------------------------------
    # free memory and stop logging
    saga_api.SG_Get_Data_Manager().Delete()
    lis_cmd.Stop_Logging(starttime)

Comparison with the LIS Pro 3D GUI Workflow

Tile-by-Tile Processing

In the Getting Started section of this tutorial series, we manually applied the processing tools to a small point cloud subset. This worked because the dataset was small enough to fit into memory.

In the Putting It into Practice section, we introduced tile-by-tile processing using the multitile-tools. We learned that processing larger point cloud datasets typically requires dividing them into smaller, manageable units with sufficient overlap to account for trees near tile boundaries.

We will follow the same approach in our Python scripts. Each script processes the point cloud tile by tile, using the processing units created during data preparation.

The scripts essentially reproduce the workflows we already explored through the LIS Pro 3D graphical user interface, including the same processing steps and parameter settings. The main difference is that we now execute these steps automatically rather than manually.

Note

All processing parameters are defined centrally in project_config.py, rather than directly in the individual processing scripts.

This means that when we want to apply the workflow to another dataset, we generally only need to adjust project_config.py, without modifying the processing scripts themselves.

Parallel Processing

Our Python scripts also support parallel processing to speed up the analysis.

Because the processing units can be handled independently, we can distribute them across multiple CPU cores instead of processing every tile sequentially.

This allows us to make better use of the available computing resources and potentially reduce processing times significantly, especially when working with larger TLS datasets.

The actual speedup depends on the number of CPU cores, available memory, and the computational requirements of each processing step.

Run All Processing Scripts Using a Runner Script

Now that we have prepared the individual processing scripts, we will create an additional Python script called run_script.py.

This script allows us to execute multiple processing steps sequentially with a single command. It also automatically creates separate log files for standard output and error messages, making it easier to monitor processing and identify potential problems.

In VS Code, create a new file named run_script.py in the scripts folder. Copy the following code into the file and save it:

from subprocess import run
import os
import sys
import argparse

LOG_DIR = "../log"


def parse_args():
    parser = argparse.ArgumentParser(description="run a script with automated logging.")

    parser.add_argument(
        "scripts",
        type=str,
        nargs="+",
        help="the script(s) to execute, e.g. 'my_script.py' or 'my_script0.py my_script1.py my_script2.py'",
    )

    return parser.parse_args()


if __name__ == "__main__":
    args = parse_args()

    if not os.path.exists(LOG_DIR):
        os.makedirs(LOG_DIR, exist_ok=True)

    for script in args.scripts:
        scriptname, fileext = os.path.splitext(script)
        std_log = os.path.join(LOG_DIR, f"{scriptname}.std.log")
        err_log = os.path.join(LOG_DIR, f"{scriptname}.err.log")

        if fileext == ".py":
            command = f"{sys.executable} -u  {script} > {std_log} 2> {err_log}"
            print("running: " + command)
            run(command, universal_newlines=True, shell=True)
        else:
            print(f"Skipping {script} (unknown file extension)")

Execute the Processing Pipeline

Now that all scripts are in place, we are ready to execute the complete TLS-based forestry processing workflow.

Open the Miniforge Prompt

First, open the Miniforge Prompt on your computer.

Activate Your Conda Environment

Activate the Conda environment containing your LIS Pro 3D installation. Replace myenv with the name of your environment.

(base) C:\Users\LiDAR>conda activate myenv

(myenv) C:\Users\LiDAR>

Navigate to the Scripts Folder

Next, navigate to the scripts folder of your project. Replace the example path with the actual location of your project.

(myenv) C:\Users\LiDAR>cd C:\LiDAR_Data\tutorials\forestry_tls\scripts

(myenv) C:\LiDAR_Data\tutorials\forestry_tls\scripts>

Run the Complete Processing Pipeline

Execute the following command to run all processing scripts sequentially:

python run_script.py 01_data_preparation.py 02_ground_classification.py 03_stem_detection.py 04_tree_segmentation.py 05_tree_models.py

The runner script will execute the individual processing steps in the specified order:

Data preparation: Create the virtual point cloud dataset and processing units.
Ground classification: Classify ground points and prepare the point clouds for forestry analysis.
Stem detection: Detect tree stems and derive their positions.
Tree segmentation: Segment the point cloud into individual trees.
Tree modeling: Reconstruct 3D tree structures and derive tree-level metrics.

The pipeline automates the complete workflow, from initial data preparation to the extraction of tree information, without requiring manual interaction with the LIS Pro 3D graphical user interface.

Note

The runner script stores standard output and error messages in separate log files in the log folder.

Since the processing scripts are executed sequentially, check these log files for potential errors before assuming that the entire pipeline has completed successfully.

Inspect the Results

Once processing is complete, navigate to the corresponding output folders in your project directory.

You can open the resulting datasets in LIS Pro 3D to inspect the classified point clouds, detected stem positions, individual tree segments, reconstructed tree models, and derived tree metrics.

Compare these results with those obtained in the earlier GUI-based tutorial sections. The processing workflow and parameter settings are essentially the same, but the individual steps are now executed automatically through Python (and for the whole plot area).

Summary

In this tutorial section, we completed the transition from the LIS Pro 3D graphical user interface to a Python-based workflow for automated TLS forestry analysis.

Building on the project configuration and data preparation steps from the previous tutorial section, we created individual Python scripts for ground classification, stem detection, individual tree segmentation, and 3D tree reconstruction. These scripts reproduce the processing steps we previously performed manually, using the same tools and parameter settings.

We also introduced tile-by-tile processing and parallelization to handle larger point cloud datasets efficiently. By distributing independent processing units across multiple CPU cores, we can make better use of available computing resources while keeping memory requirements manageable.

Finally, we created a runner script that executes the complete processing pipeline in the correct order and automatically generates log files. This allows us to run the entire workflow, from data preparation to the extraction of tree-level information, with a single command.

With this approach, we have established a reproducible, scalable, and largely automated workflow for TLS-based forestry analysis. By keeping processing parameters in a central configuration file and separating individual processing steps into dedicated scripts, the workflow can be adapted to new datasets with relatively little effort.

Together, the tutorials in this series demonstrate how LIS Pro 3D can be used for both interactive exploration of individual trees and automated analysis of larger forest areas.

< Previous section Next section >